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the distribution P[h(m)] of functions of magnetizations. The order parameter is interpreted as the 
histogram of probability distributions of individual magnetizations. In the limit of zero temperature 
and small transverse fields, to leading order in F magnetizations m ~ become relevant in addition 
to purely classical values of m ps ±1. Self-consistency equations for the order parameter are solved 
numerically using Quasi Monte Carlo method for K = 3. It is shown that for an arbitrarily small 
T quantum fluctuations destroy the phase transition present in the classical limit F = 0, replacing 
it with a smooth crossover transition. The implications of this result with respect to the expected 
performance of quantum optimization algorithms via adiabatic evolution are discussed. The replica- 
symmetric solution of the classical random Tf-Satisfiability problem is briefly revisited. It is shown 
that the phase transition at T = predicted by the replica-symmetric theory is of continuous type 
with atypical critical exponents. 
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I, INTRODUCTION 

Quantum phase transition (QPT) is a transition be- 
tween different ground states driven by quantum fluctua- 
tions and controlled by certain parameters, for example, 
an external magnetic field. Study of QPTs in systems 
with strongly interacting spins attracted attention in the 
field of quantum computing due to the possibility of cre- 
ating massively entangled states at the quantum criti- 
cal point [l| and the relevance of QPTs to the analysis 
of the performance of quantum algorithms for solving 
classical combinatorial optimization problems (COPs) 
0, H, 0, H, @|. Quantum mechanics offers an alterna- 
tive to the mechanism of thermal fluctuations for the 
transitions between the states, which can be exploited 
in the optimization procedures 0, @|. QPT in this pa- 
per will be studied in the context of a general-purpose 
quantum adiabatic algorithm (QAA) proposed by Farhi 
and coworkers @|. In its simplest form the algorithm is 
defined via a quantum TV-spin Hamiltonian that is a sum 
of two terms 

N 

H = n cl (a^...,a z N )-rJ2^- (1) 

i=l 

The first operator term is derived from a cost (energy) 
function of classical spins H c i(si, . . . , sjv) by replacing 
each classical spin Si — ±1 with a Pauli matrix, of. The 
ground state of this operator encodes the solution of a 
classical COP described by TL C \. The second term de- 
scribes spin coupling to the external magnetic field oc T 
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applied in the transverse direction (e.g. along the posi- 
tive x axis). At the start of the algorithm, T is made very 
large and the ground state of H (0) is prepared with all 
the spins pointing in x direction. Then T=r(i) is slowly 
reduced to zero while the state of the quantum system 
remains close to the instantaneous adiabatic ground state 
of H(t) — provided that the condition (* |^-ff|*o) < 
{Ei-Eo) 2 is satisfied. Here H(t)\^ n (t)) = E n (t)\^ n (t)) . 
At the end of the algorithm at T = 0, the system is found 
in a state which is a superposition of spin configurations 
corresponding to all degenerate global minima of H c \. 
The runtime of the algorithm is proportional to l/g 2 lhl , 
where g m i n = minr(£'i —Eo) is a minimum of the energy 
gap [9] taken over the range of T. 

It has been noticed several decades ago that properties 
of the solution space of complex COPs are closely related 
to those of spin glass systems ficl . It has been also 
recognized |12| that many of the spin glass models are in 
almost one-to-one correspondence with computationally 
hard COPs encountered in practice and forming a class 
of NP-hard problems. 

Whereas theoretical computer science is mostly con- 
cerned with the worst-case complexity, from the statisti- 
cal physics perspective the main interest lies in the typi- 
cal running time of algorithms over the random ensemble 
of p roblem instances (or samples of spin glass system) 
[ill . fl4 |. When this expected runtime scales exponen- 
tially with the number of spins, the COP is considered 
intractable. This intractability was linked to so-called 
threshold phenomena flR [l6l . [ItJ in NP-complete prob- 
lems. In physics community, these threshold phenomena 
were recognized as phase transitions in models of classical 
spin glasses [l8j]. Many NP-complete problems, includ- 
ing the most basic of them — random if -Satisfiability (or 
if-SAT) — correspond to infinite-range dilute spin glass 
models with if-local interactions, i.e. H. c i(si, . . . , sjy) is 
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given by a sum of interaction terms; each involving a 
set of K spins, chosen at random from a set of size N. 
In contrast to finite-dimensional models, the topology of 
links corresponding to spin couplings is completely ran- 
dom, with no non-trivial correlations. The random en- 
sembles of instances are described by a single parameter 
the connectivity 7 which is the number of interaction 
terms per spin, 7 = M/N . The probability for a given 
spin to be involved in d interactions is Poisson with the 
finite mean value of d equal to Kj. This is different 
from infinite-range fully-connected spin models such as 
the Sherrigton-Kirkpatrick model [19(|, where the value 
of d = N — 1 scales with the number of spins. 
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FIG. 1: Thick black (1) and gray (2) lines show two possi- 
ble forms of quantum phase diagrams on transverse field V 
vs connectivity 7 plane for random K-SAT problem. Black 
line (1) corresponds to the quantum dilute ferromagnet. Gray 
filled rectangle shows the region of interest in this paper with 
small transverse fields, F « 1. Dot-dashed line depicts the 
scaled exponent a = a(~f) of the median runtime T of a clas- 
sical algorithm, T ~ exp(— aN), over an ensemble of problem 
instances with the same 7. 

Classical infinite-range spin glass models in the dilute 
limit have been studied, though recent results concen- 
trate on the zero-temperature limit [10, HH, H3, El, 
[25l . HE]- However, very little is known about their re- 
spective quantum versions [described by the Hamilto- 
nian ([I]) with H c \ corresponding to an infinite-range di- 
lute spin glass] despite a lot of interest in these models 
from the perspective of quantum computing. Recently, 
quantum versions of random Exact Cover and other re- 
lated optimization problems have been studied us- 
ing a generalized annealing approximation [1H|- At the 
same time, fully connected infinite-range quantum spin 
models have been analyzed in the literature using vari- 
ous approximations. This includes quantum versions of 
the Sherringtori-Kirkpatrick [H, H3, [HI, E3, random 
Heisenberg [3] , p-spin and random energy models [35| . 
Exact solutions in quantum spin glasses are mainly lim- 
ited to one-dimensional models [3a. [37I. 



Numerical studies [16] have demonstrated that the typ- 
ical runtime of known classical algorithms applied to en- 
sembles of randomly generated instances of AT-SAT and 
similar models, as a function of 7, peaks at the point 
of static transition, which is a major bottleneck of clas- 
sical optimization algorithms (see Fig. [J). This can be 
understood by analogy with the critical slowing down 
of the dynamics in the vicinity of phase transitions in 
problems without disorder. Similarly, we expect that the 
dynamics of QAA for random if-SAT could be governed 
by the corresponing quantum phase tranisiton (QPT). If 
the system underwent a QPT as the value of T is lowered 
from a large value to 0, the gap would attain its minimum 
value at the point of the transition. The critical exponent 
associated with the singularity of the free energy would 
determine the scaling of the minimum gap (which would 
have the form of an exponential or stretched exponential 
0). 

In this paper we concentrate on static transition that 
corresponds to satisfiability transition at zero tempera- 
ture. We concentrate on K = 3 as the most interesting 
case. It is the smallest value of K for which if -SAT 
is NP-complete. Moreover, random if-SAT undergoes 
random first-order phase transition for K ^ 3. As it is 
the case with all random first-order transitions, the static 
transition is preceded by the dynamic transition. Results 
for similar model — if-XOR-SAT, or dilute p-spin glass 
— at finite temperature indicate [H] that the free en- 
ergy remains analytic across dynamic transition, which 
would imply that the static transition is the real bottle- 
neck of simulated annelaing algorithm. While giving cre- 
dence to the idea of the analysis of static transition, this 
picture may not necessarily apply to fC-SAT for K — 3, 
where the dynamic transition is accompanied by another, 
condensation, transition [H, Hfil- Due to difficulties of 
replica-symmetry-breaking analysis in quantum case, we 
have only performed the replica-symmetric analysis. Al- 
though replica-symmetric approximation is capable of 
correctly capturing the existence and qualitative prop- 
erties of static transition, it fails to describe the dynamic 
transition and overestimates the critical threshold 7 C . 

In Fig. Q] we sketch two conjectured forms of QPT line 
r = r c (7). One possiblity r c (7) changes continuously 
from the value of at 7 = j c . Alternatively, it may ex- 
hibit a finite jump (i.e. r c (7 c ) = r c0 > 0) as in dilute 
transverse Ising models without frustration [39l.[4fj|. An- 
other (third) possibility is that the phase transition at 
T = disappears for any finite T > 0. One may distin- 
guish between these cases by setting T <Cl and studying 
the free energy for a range of values of 7 containing j c , 
as shown in Fig. [TJ In QAA the parameter T(t) decreases 
with time, corresponding to a vertical line in the (7, T) 
plane as shown in Fig. [H The central result of this paper 
is that it is the third possibility that takes place: quan- 
tum effects (the transverse field T) in the QAA Hamilto- 
nians (U) make the static phase transition disappear; the 
free energy becomes analytical in the vicinity of 7 C for 
small but finite T. 
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It should be mentioned in passing that certain highly 
symmetric examples of COPs have been constructed 
[ill, where a total spin is an exact quantum number 
of the Hamiltonian H of Eq. ([!]) and QAA fails due to 
the onset of a large spin tunneling through a broad, order 
N , semiclassical barrier with amplitude that scales down 
exponentially with N [H, Hj3]. However, in spin glasses, 
quantum evolution does not correspond to large spin dy- 
namics. Instead, an exponentially large (in N) number of 
deep local minima of the classical energy are connected 
by an extremely large number of tunneling paths with 
amplitudes proportional to high powers of T. This pic- 
ture as well as the analysis of QPTs is more relevant for 
understanding the typical complexity of QAA for NP- 
hard problems such as if -SAT. 

This paper is organized as follows. Section Ull presents 
a brief overview of important results for the classical ver- 
sion of random if-SAT and discusses the relationship be- 
tween the present work's replica-symmetric analysis of 
quantum if -SAT and that of the classical if-SAT cor- 
responding to the limit of T = 0. We formulate the 
quantum version of if -SAT and analyze it using replica- 
symmetric theory in Sec. lIIIl This is followed by the anal- 
ysis of small magnetic fields T in Sec. IIVI In Sec. [V] we 
revisit the classical T = random if -SAT to demonstrate 
that the replica-symmetric analysis predicts a continuous 
phase transition; it was previously thought to be of ran- 
dom first-order type. In Sec. IVII we present the numer- 
ical results for both finite-temperature classical if -SAT 
and zero-temperature quantum if-SAT. We concentrate 
on if = 3, which is the most intersting case. Since we 
utilize the replica-symmetric approximation in the anal- 
ysis of quantum if-SAT, we compare these results with 
those predicted by the replica-symmetric theory for fi- 
nite temperature classical if-SAT (despite the fact tools 
to study replica symmetry breaking in classical if-SAT 
have appeared recently). In the Conclusion we discuss 
our results, especially in relation to quantum adiabatic 
algorithm and describe possible extensions of the present 
work. The mathematical details of the calculation of 
the replica free energy functional are relegated to Ap- 
pendix[A] Appendix[B]discusses correspondence between 
the replica-symmetric ansatz and the Bethe-Peierls ap- 
proximation. A novel Quasi Monte Carlo algorithm used 
in numerical calculations is described in Appendix ICl 

II. CLASSICAL STATISTICAL MECHANICS OF 
RANDOM if-SAT: MONASSON-ZECCHINA 
REPLICA SYMMETRIC SOLUTION AND ITS 
CONNECTION TO THE PRESENT WORK. 

An instance of random if-SAT is a system of N clas- 
sical spins with the energy function that is written as a 
sum of M terms: 

W c i(si, sn) = E(s il ,...,Si K ;J e ). (2) 

e=(ii...i K )<£ 



Each term is associated with a if -tuple e = . . . , ix)- 
If spins labeled by i — 1,...,N are viewed as vertices 
of some graph, if -tuples e correspond to its hyperedges. 
The set of all hyperedges for a given instance is labeled 
£. Hyperedges corresponding to each term are chosen 
independently and uniformly at random; hence with each 
instance of random if-SAT we may associate a realization 
of a random hypergraph. This represent the geometric 
part of disorder. 

Each term defines a constraint involving spin variables 
s% x , ■ ■ . , Si K . The cost function W c i(si, . . . , sk) can be 
either zero or some positive value representing the energy 
penalty for those combinations (s%, . . . , sk) that violate 
the constraint. 

For if-SAT the constraints penalize exactly one out 2 K 
assignments. The cost function is chosen in the following 
form 

E J (s u ...,s K ) = 2f[ 1 -±^. (3) 

Here J = (Ji, . . . , Jr-), where Jg = ±1, denotes the 
combination of if spin values that is assigned an energy 
penalty of 2 [65| . The argument J of the cost function 
will be written as a subscript unless it refers to a specific 
hyperedge as in Eq. ((2]). The values of disorder variables 
J e are chosen independently and uniformly at random for 
each constraint. The corresponding probability distribu- 
tion assigns the probability of 1/2 K to each realization 
of J: 

pw=f ? (j '-":* (j ' + " - w 

The energy |f5]) equals twice the number of violated 
constraints. When the number of constraints M is suf- 
ficiently small, all of them may be satisfied at the same 
time and the energy is zero. The properties of random 
if-SAT are studied in the limit when the number of vari- 
ables N and constraints M goes to infinity, while the 
constraint-to- variable ratio 7 = M/N is kept constant. 
In this limit the fraction of variables involved in d con- 
straints is Poisson with mean if 7 

f4Kj) = ^(if7) d e-^, (5) 

so that each variable appears in if 7 constraints on aver- 
age. 

It has been shown by computer studies that there exists 
a threshold j c such that with overwhelming probability, 
there exists a configuration of N spins with zero energy 
if and only if 7 < j c (in the limit of large N) . In the lan- 
guage of statistical mechanics, the random if-SAT un- 
dergoes a phase transition between the satisfiable (SAT) 
and unsatisfiable (UNSAT) phases at 7 ~ j c . The inter- 
action term (J3j) imposes a "weak" constraint on the spins 
involved in it. For this reason, unlike the Viana-Bray 
model with Ising interactions, the phase transition for 
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random K-SAT does not coincide with the percolation 
transition for the corresponding hypergraph. For 3-SAT, 
the percolation transition takes place at 7 pC rc = 1/6, 
while the "experimental" value of the satisfiability thresh- 
old is 7 C w 4.2 fl6l |. The exact value of 7 C for random 
K-SAT for K > 3 is not known. 

Random if-SAT can be formulated as a statistical me- 
chanics problem by introducing the artificial temperature 
T = 1/(3 and writing the Gibbs free energy 
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(6) 



The extra factor of 1/N ensures that this is the free en- 
ergy per spin so that F does not scale with N. It is 
related to the total internal energy E and the total en- 
tropy £ via the standard identity: 



F= — (E-TE) 



(7) 



In the limit T = thermal fluctuations disappear and the 
second term in Eq. J7]) vanishes. In this limit E converges 
to the minimum value of energy H c \. Therefore, F = for 
7 < 7c. Note that in random K-SAT there is no region 
where the minimum number of violated constraints is 
o(N) except in the immediate vicinity of j c . For 7 > 7 C 
this number is O(N) and F > 0. 

Instance-to-instance fluctuations of F are small: o(l). 
Therefore, with overwhelming probability a randomly 
chosen instance has free energy within o(l) from (F), 
which is the disorder-averaged value. This is the central 
quantity which is computed using the replica method. 
We briefly discuss the main results obtained in [10, HH- 
The authors demonstrated that the disorder-averaged 
free energy of random if-SAT corresponds to the ex- 
tremal value of the free energy functional 

T[P(h)] = 

/OO (>OG 
d/ii-- - / dh K P(hi)...P{h K ){Uj(h 1 ,...,h K ))j 
-OO J — OO 

Here (■ ■ - }j denotes averaging over the parameters Ji — 
±1 (I — 1, . . • , K) with equal weights assigned to all 2 K 
possibilities. The function Uj({he}) is defined as 

Uj(h\, . . . , h K ) = 2 min (l, {J\h 1 ) + , . . . , (J K h K ) + ) . 

(9) 

Here and throughout the paper we use a shorthand (■■■)+ 
which we define as follows 



(*)+ 



x for x > 0, 
for x sC 0. 



(10) 



The function P(u>) in ([8]) is the Fourier transform of the 
distribution P(h): 



dhe~ iwh P(h). 



(11) 



The function P*{h) is found by extremizing 6J-[P(h)] 
subject to the constraint f_°° dhP(h) — 1 a has the 
meaning of the histogram of effective fields hi associated 
with each spin. Whenever hi 7^ spin Si takes the same 
value Si — sgn hi in all spin configurations with the lowest 
energy. The absolute value \hi\ is one-half of the energy 
cost needed to flip it. 

The fraction of frozen (s.t. hi 7^ 0) spins q = 
J ^ dhP*(h) + J^g 00 dhP*(h) is the order parameter as- 
sociated with the satisfiability transition. In the sat- 
isfiable phase q — 0, corresponding to P*(h) = 6(h), 
whereas the unsatisfiable phase is described by finite 
q > 0. 

The simplest solution P*(h) of the extremality condi- 
tion for the functional J51) is 

+00 

P(h)= J2 e- K ^ K - 1 I w (K 7 (q/2) K -l)6(h-k), 

k— — oo 



where Ik(x) is the modified Bessel function of first kind. 
The value of q may be determined self-consistently from 



h)of 

m 
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(12a) 



For K = 3 and 7 > j d ~ 4.667 Eq. (fl2a|) has two stable 
solutions: the trivial q = and the non-trivial q > 0. 
The non-trivial solution does not becomes stable until 
7 > 7c w 5.181. The corresponding bound is very close 
to the annealed bound of 7 an n = ln2/ln(7/8) « 5.191 
and greatly overestimates the "experimental" value of the 
satisfiability threshold 7 cxp ~ 4.2 from computer simula- 
tions [16]. 

A similar integer-delta-peaks solution [44j for the or- 
der parameter in the Viana-Bray model [45J| was shown 
to be unstable in the longitudinal sector (i.e. within the 
replica-symmetric ansatz)[46j]. The longitudinally stable 
solution exhibited a continuous part in addition to delta- 
peaks. Though the appearance of the continuous compo- 
nent is believed to signal the breakdown of replica sym- 
metry, the replica-symmetric result may still be useful if 
regarded as a type of variational approximation. 

The incorporation of the continuous component led to 
an improved upper bound of the satisfiability transition 
7c w 4.60 obtained numerically [2ll |. This problem will 
be revisited in Sec. |V] and we will show that although 
the value of 7 C had been determined correctly, the phase 
transition predicted by the replica-symmetric theory is 
actually continuous rather than first-order as was claimed 
in Ref. [HI]. 

Subsequent analysis by M. Mezard and R. Zecchina 
of 1-step replica symmetry breaking (RSB) in random 
if-SAT improved the bound for satisfiability threshold 
to 7c w 4.267 [H,[23|]. It is believed that this 1-step RSB 
solution is stable. What made the T = RSB analy- 
sis tractable (and yet required a lot of numerical effort) 
was the integer-delta-peaks ansatz for the distribution 
of effective fields within each pure state. It is a daunting 
task to extend 1-step RSB analysis to finite temperatures 
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(where non- integer effective fields are certain to exist), let 
alone including quantum effects. This paper only consid- 
ers the replica-symmetric solution. 

Using replica-symmetric analysis to study the quantum 
problem may have some merit. It has been argued in the 
literature UtJ, based in part on results on quantum SK 
model that effects of quantum tunneling may sta- 

bilize the replica symmetric solution. Even if true, such 
symmetry must break down for extremely small trans- 
verse fields r = o(N)/N or in the limit T/T <C 1. Indeed, 
the purely classical limit T = should be described by 
the 1-step RSB solution obtained in Ref. [23] . 

III. REPLICA SOLUTION OF QUANTUM 
A"- SAT 

A. Replica-symmetric free energy functional 

The quantum Hamiltonian given by Eq. (jTJ is a sum 
of two terms: the purely classical term describing the in- 
teraction of Ising spins and the quantum term describing 
the coupling to the external magnetic field applied in the 
transverse direction. By employing Suzuki- Trotter trans- 
formation, the problem of finding the partition function 
Z = Trc~P H can be reformulated as that of computing 
the partition function of the purely classical model. The 
corresponding classical partition function is written as a 
sum over all possible paths Si(t): 

[{•♦(«)}] 

where the functional JC[s(t)] is given by 

£[»(*)] = — ln(tanhTAt) ^ s(t)s(t + At) 

t=Q,At,...,/3-At 

+ ^ln(-sdnh2rAt). (13a) 

The time variable t takes L discrete values t = kAt (At = 
0/L). Periodic boundary conditions Si((3) = s,(0) are 
assumed. 

The sum (|13p is over N x L spin variables labeled by 
i = 1, . . . , N and t = kAt. In anticipation of the limit 
L — > oo that will be taken eventually, we treat time as a 
continuous variable. In particular, we write ji 3 dt ■ ■ ■ to 

mean J2k=o ^ ' " • We use square brackets for writing 
functionals and for indicating sets labeled by continuous 
variables. Sets indexed by a discrete variable will be des- 
ignated using curly braces. To avoid ambiguities we may 
adorn brackets or braces with subscripts and superscripts 
to indicate index variables and ranges (e.g. [{si(t)}i]f_ ) 
A constant in expression (|13a[) ensures proper normal- 
ization of the statistical sum (jT3j) . It can be verified 
that CEU) reduces to Z^> = (2cosh/3r) A ' for the non- 
interacting problem {H c \ = 0). 



We choose to write the classical Hamiltonian ([2]) in the 
following form 

"^cl ^ ^ c ii ...IK ^ 5 ■ ■ • J , Jii ...Ik )j (14) 

U<«2<-<i/f 

where the cost function Ej(s\, . . . ,Sk) for K-SAT is 
given by Eq. J3]). Disorder variables Ji 1 ...i K are assumed 
to be uniformly distributed according to Eq. (|4]). 

The value of (H 1 ...i K is chosen to be 1 if the instance con- 
tains a constraint involving a set of variables i\,...,iK, 
and zero otherwise. Random variables Ci 1 ...i K are statis- 
tically independent and distributed according to 

P^ = {}-^k)^) + ^k^-^ (15) 

In the asymptotic limit (N — » oo) the number of con- 
straints will be M = 7./V. Form lfT4|) is preferable to ([2]) 
because it emphasizes the long-range character of random 
K-SAT. 

In this paper we will keep the derivation as general 
as possible. Formulae written without expanding J3]) 
will be — by substituting appropriate expressions for 
Ej(s\, . . . , sk) and p(J) — directly generalizable to any 
random combinatorial optimization problem with binary 
variables and if-local interaction (e.g. iT-XOR-SAT, 
if-NAE-SAT, 1-in-A- SAT). 

The central physical quantity of interest is the disorder 
averaged value of the free energy (F) — — (In Z) . This 
is the same as the value of the free energy for a typical 
realization of disorder, the free energy (in contrast to 
Z) being a self-averaging quantity. We use the replica 
method to perform the disorder averaging. The average 
of the logarithm is rewritten using the following identity: 

(hxZ) = lim ^-(Z n ). (16) 

For integer n, Z n is the partition function of a system 
of n non-interacting replicas of the original random in- 
stance. Computing (F) will require performing the ana- 
lytical continuation in n. The gist of the method is that 
disorder averaging in the expression for (Z n ) is done prior 
to performing the sum over classical spin configurations. 

[{»?(*)}] 

.(17) 

where the replica index a runs from 1 to n, effectively 
increasing the number of spin variables to N x L x n. 

Disorder averaging couples together formerly non- 
interacting replicas. However, it also transforms the di- 
lute model with strong 0(1) interactions into a com- 
pletely connected model with weak 0(1/N K ~ 1 ) interac- 
tion. This permits the exact evaluation of the sum over 
the spin variables using mean field theory. We express 
the mean field solution in terms of a set of order param- 
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eters: spin correlation functions 



In the thermodynamic limit, the partition function pT 
can be written in the form of a functional integral: 



(18) 



(Z n ) = J VQVXe 



■Nnf3T[{Q},{\}] 



(19) 



The argument of the exponential is (up to a factor) 
the free energy functional T that depends on correla- 
tion functions {Q ai ...a p (ti, ■ ■ ■ ,t p )} as well as Lagrange 
multipliers {A ai ... 0p (ii, • • • ,t p )} that enforce constraints 
((18)1 . In Eq. ((19)) we have suppressed indices and time 
arguments for conciseness; similarly VQ and PA are a 
shorthand for multiple functional integrals. 

In the limit N — > oo the integral lfl9)) is dominated by 
the saddle-point value of T: 



F = — 



1 



iVn/3 



ln(Z")=P[{Q*},{A*}] 



(20) 



The right hand side is evaluated for {Q* a {t)}, {A„(£)} 
that make T stationary with respect to small variations. 
Note that in the following we will use a calligraphic T to 
indicate a functional and a roman F to denote its value 
at the saddle point. 

In practice, working with an infinite set of time- 
dependent correlation functions is infeasible. Instead, 
as often done in the analysis of quantum spin glasses 
[30l l34j. we resort to the static approximation. We solve 
stationarity condition for the reduced set of functions 
— those that are independent of time arguments. Note 
that consistency requires that if Q ai ...a p (ti, . . . ,t p ) are 
replaced by their static counterparts Q ai ...a p > any time- 
dependence be ignored for X ai ...a p (ti, ■ ■ ■ ,t p ) as well. Im- 
plemented in this form, the static approximation may be 
regarded as a type of variational approximation. 

Integrating out {A 0l ... }, we may write T{{Q}) as 
a function of {Q ax —a p } alone. It may be verified that 
static Q ai ...a p are the time-averaged dynamic correlation 
functions: 

Qax...a p = -Tp J dti . . ■dtpQ ai ... ap (ti, . . . ,t p ). (21) 

We work within replica-symmetric ansatz, which posits 
that Q ai ...a p at the saddle-point of T are symmetric with 
respect to permutations of replicas. Due to this symme- 
try, not all Q ai ...a p are independent. The value oiQ ai ,„ a 
may only depend on the set of numbers fei, fe, • • • which, 
respectively, indicate the number of distinct replica in- 
dices that appear exactly once, twice, etc. We will write 



Qa 1 a 2 ...a kl fci fci&2&2 •■ -&fc 2 bfc 2 ■•• — Qfeifc 2 .. 



(22) 



where {a,}, . . . are all distinct. Although for finite 
integer n, the inequality ^ r k r ^ n must hold, perform- 
ing the analytical continuation to n — > requires the 



knowledge of (Z n ) for all integer values of n. Thus, para- 
doxically, in the limit n — > 0, the values k r may run from 
1 to 00. 

Note that in classical limit T — 0, only two paths 
[s(t) ee +1 and s(t) ee -1] contribute to fI7)l. Due 
to that, the static approximation becomes exact in this 
limit, and the order parameters Q{k r } may depend only 
on p = Y^r^r+i as evidenced from Eq. (jTSj) . It has 
been recognized in the analysis of classical Viana-Bray 
model by I. Kanter and H. Sompolinsky that the or- 
der parameters Q p are the moments of the probability 
distribution P(m) of average spin magnetizations. For a 
quantum model, Q{k r } are related to the functional dis- 
tribution P[h(m)], where functions h(m) are defined on 
the interval [—1; 1]: 



Q {kr} = \ [dh(m)} P[h(m)} [] 



J dm* 



-0h(m) r \ k 



Id 



m e 



-l3h(m) 



(23) 

That the right hand side of (|23)l is a functional integral 
is indicated by the use of square brackets (J[dh(m)] • • •). 
Such notation is customary in quantum field theory (see 
e.g. [HI]) and is consistent with our practice of using 
square brackets to indicate sets indexed by continuous 
variables. Regular multidimensional integrals will be 
written using curly braces (e.g. J {dmi}^ =1 ■ ■ ■) . Note 
that integrals over magnetizations run from —1 to +1. 

We refer to functions h(m) as effective fields. It can 
be guessed from the formm of f23)) that these effective 
fields represent probability distributions of individual 
spin magnetizations viapi(m) oc e~P h ( m \ The distribu- 
tion P[h(m)] is the histogram of effective fields hi(m) as- 
sociated with each spin. It may be interpreted as a prob- 
ability distribution of probability distributions of magne- 
tizations. Such constructs appear in replica analysis of 
classical problems in the description of replica symmetry 
breaking (RSB). As one can see, in the quantum case 
they are already present at the replica-symmetric level. 
Note that the effective fields h(m) are defined only up to 
a shift by an arbitrary constant h(m) — > h(m) + const. 

We express F({Q{k r }}) m terms of the distribution 
P[h(m)] as a sum of two terms, which we will call a 
"quasipotential" V and a "quasientropy" S; themselves 
dependent on P[h(m)]: 



nP[h(m)]j = 7 V[P[Mm)]] - S\P[h{m) 



(24) 



We have used double square brackets to indicate that ar- 
guments of J 7 , V and S are functionals. Detailed deriva- 
tions are given in Appendix [Aj here we provide the re- 
sulting expressions. For the quasipotential V[P[/i(m)]] 
we obtain 

K 

V= [{d^(m)}jy Y[P[h e (m)} x (Uj[{h e (m)}})j, 
J 1=1 

(25) 

where (• indicates averaging over 2 K possible re- 
alizations of vector J. The functional integral over 
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hi(m), . . . , hx (m) describes averaging over probabil- 
ity distributions P[hi{m)] of the quasipotential density 
Uj[h\{m), . . . , hxim)] given by the following expression: 



1 K r 
Kj[{ht(m)}] = -^ln/ dme 

" 1=1 J 



In / {dmi} e=1 e 



K -pEj(m 1 ,...,m K )-f3Y.Y = ihi(m l ) 



(26) 



Integrals over magnetizations run from —1 to +1. We 
write / {dmi}jf =1 • • • to indicate the if-dimensional inte- 
gral over magnetizations mi, . . . , mx- 

The function Ej(m\, . . . , tuk) that appears in Eq. lf26|) 
is multilinear in mi, . . . , mjf and coincides with 2?j(. . . ) 
when {mi} G {±1} . These two conditions determine 
it uniquely. For K-SAT the expression is obtained by 
formally replacing discrete spin variables in Eq. ([3]) with 
continuous magnetizations {mi}: 



Ejimi, . . .,m K ) 



1 + Jrmi 1 + J K m K 



It is easily seen that for any t one may write 
Ej(mx, . . . ,m#) = Ag + Bimi, where Ai and Bi are 
independent of mi but depend on J and other magneti- 
zations {mi'}i'^i. 

For the quasientropy 5[P[/i(m)]], we obtain the fol- 
lowing expression: 

S= J [Ah{m)]C[h{m)\ J [dw(m)] dm ^ m)Hm) S[uj(m)}, 

{28]) 

with S[co(m)] given by 
£ = P[u>(m)] ( 1 — i / dmu(m)«o(m) — In P[u(m)] I , 



(28a) 

which in turn is written in terms of the functional 
Fourier transform of P[h(m)} that we denote P[u)(m)}. 
It is implied that the normalization inside the func- 
tional integral over cj(m) is such that the inverse 
Fourier transform of P[uj(m)] reproduces P[h(m)], i.e. 
/[dcj(m)]e i / dmw ( m ) /l (" l ) J P[w(TO)] = P[h(m)]. 

The functional C[h{m)} is given by the following ex- 
pression: 



C[h(m)] 



dme- f3h(m \ 



(29) 



The function uo(m) that appears in Eq. (|28ap is en- 
tirely due to the kinetic term /C[s(t)]. In the limit of 
continuous magnetizations {L — > oo) it can be evaluated 
in closed form: 



-0ua(m) _ 



(3T 



\/l — m 2 



h (pry/ 1 - m 2 ) 

+ 5(m-l) + 5(m + l). (30) 



Observe that in the limit r = only contributions from 
m = ±1 are expected. We demonstrate in Appendix I A 41 
that the free energy functional ([24| may be re-expressed, 
using the reduced order parameter P(h), in the form 
given by Eq. JSJ. 

It would seem from the form of Eq. lfT9|) that the free 
energy should correspond to the minimum of the free en- 
ergy functional ([24|l . Because of the peculiar nature of 
the limit n — > 0, this is not the case. In Appendix I A 41 we 
show that in the classical limit (r = 0) the free energy 
is a local maximum with respect to symmetric pertur- 
bations of P(h) [i.e., such that 8P(~h) = SP(h)] and 
a local minimum with respect to antisymmetric pertur- 
bations [s.t. 5P(—h) = —SP(h)]. The quantum case is 
considerably more complex; fortunately, we only need to 
make sure that P[h(m)] makes the free energy functional 
T stationary and do not care whether it is a minimum 
or a maximum. 

A few notes must be made about approximations made 
in this section. The assumption of replica symmetry is 
justified for sufficiently small connectivities 7; above the 
replica-symmetry-breaking transition (7 > 7rsb), it be- 
comes an approximation. In contrast, the static approx- 
imation is not guaranteed to be exact anywhere except 
r = 0. It is a type of mean-field approximation, whereby 
fluctuating spins are replaced by average magnetizations. 

The physical interpretation of the static approximation 
is rather intuitive. One can define the effective classical 
model with discrete spins replaced by continuous mag- 
netizations mi e [— 1; 1], For a specific realization of 
disorder 



Z({J})= /{dm^ilie 



-0n B!{ ({ mi };{J}) 



where the effective Hamiltonian H c s({mi}; { J}) is 
H e ff = Ej(m h ,. . .,m iK ;J il ...i K ) + ^u (mi), 

^ (31a) 
where i K )' ' ' denotes sum over all hyperedges 

Cix...i K — 1- Magnetizations m« roughly correspond to 
expectation values (of), Eq. IpTj) depends on T indirectly 
through form of Uo(m). 

In Appendix [B] we demonstrate that the replica- 
symmetric static solution is equivalent to the Bethe- 
Peierls approximation jH3] of the effective classical model 
defined by Eqs. {0TJ). 



B. 



Stationarity condition and the Monte Carlo 
method 



To complete the derivation of the replica free energy 
we need to find P[h(m)] that makes the free energy func- 
tional JF[[P[/j(m)]] stationary with respect to small varia- 
tions; its value will be the desired free energy F, formally 
a function of 0, T, and 7. The stationarity condition may 
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be written as follows: 

ST 5V 

1 



ss 



SP[h(m)] 1 SP[h{m)} SP[h{m)} 



const . 



(32) 



The arbitrary constant appears on the right hand side of 
Eq. (|32|) is a Lagrange multiplier associated with the nor- 
malization condition J dh(m) P[h(m)] — 1. Substituting 
expressions (|25| and lj28|) we will formulate the equa- 
tion that must be satisfied by the saddle-point value of 
P[h(m)]. Due to a remarkable cancelation we will able to 
write this self-consistency equation in a relatively simple 
form. 

Due to the specific form of the functionals (|25)l and 
([26|) , we may express the variation of V in the following 
form: 



SV 



SP[h{m)] 

K\ I [du(m)] Q[u(m)]£[h(m) + u(m)] - C[h{m)} 



(33) 



This identity can be used as a definition of a new func- 
tional Q[u(m)}. It is necessarily normalized to unity 
(J[du(m)]Q[u(m)] — 1). We will see that its mean- 
ing is that of the probability distribution of u{m) — 
uj (to; [{he(m)}^ =2 \) where 



uj(m; [{ht(m)}]) = ^ ^X! ln J dme 



-pht{m) 



In / {dm e } e=2 e 



K -/3Ej(m,m 2 ,...,m K )-l3Y2e=2 fl e( m e) 



(34) 



and under the assumption that J is uniformly distributed 
and /i2(m), . . . , hx (to) are taken from P[h(m)]. 

On the other hand, the variation of the quasientropy 
with respect to P[h(m)) reads 



SS 



5P[h{m)] 



i J [du(m)] C[h{m) + u(m)\ J [dw(m)] e^ raw(m >"W 
x I lnP[oj(m)] +i / dmui(m)uo(m) I . (35) 



Combining Eqs. lj33|) and (|35| uncovers the following 
system of self-consistency equations: 



K 

Q[u(m)}= / [{d^(TO)}f =2 ] Y[P[h t (m)] 

J 1=2 

x {S[u(m) - uj(m; [{h e (m)}])]) j , (36a) 



P[h(m)} = / [duj(m)]e^ dm ^ m ^ h(m '>' U0 ^ 



;d«(m)] e^ dmuj{m)u{m) Q[u(m)} 



x exp — 1 



(36b) 

In (|36a|) we use a functional generalization of 
the delta function, defined so that F[x(m)] — 
J[dy(m)] F[y(m)]S[x(m) ~ y(m)]. Note that Eq. (|36b|> 
may be written in an alternative form by expanding the 
exponential in the integrand (the term corresponding to 
d = is e- K ~<5[h{m) - u (m)j): 



P[h(m)] = J2f d (K 1 ) [{d« fc (m)}Li] Il0[«*(m)] 



d=0 



X 5 



fc=l 



ft (?7l) - U (to) - ^2 U k( r 



k=l 



pBT i 



The appearance of the Poisson distribution fd(a) = 

^i-e _a is intimately related to the hypergraph model that 
we study, as it is the distribution of the degrees (number 
of incident h ypere dges) of the vertices. From the form of 
Eqs. (|36ap . l|36bip it is apparent that h(m) are properly 
associated with the vertices of the random hypergraph, 
whereas u{m) correspond to its hyperedges. This link is 
explained in Appendix [Bj 

The system of equations |36|) can be solved itera- 
tively. Starting from some initial distribution P^ [h(m)], 
we may compute a sequence of {Q( r '[u(m)]} and 
{P (r) [ft(m)]} by applying (|36a|) and p6b]) . The limit- 
ing distribution 



P*[h{m)] = lim P {r) [h{m)] 



(37) 



must be a solution to the stationarity condition (|32|) . The 
value of the free energy is obtained from F = jV — S, 
where the quasipotential V is found by substituting 
P*[h(m)] into l(25|) . and the expression for the quasientr- 
topy S is rewritten using self-consistency equations ([36|l : 



S = K~f I [dft(m)dti(m)] P*[h(m)]Q*[u(m)] 
x (£[h{m)] - C[h(m) + u(m)]) 

[dft(m)]P*[ft(m)]£[ft(m)]. (38) 



The iterative procedure described above lends itself to 
a straightforward implementation using a Monte Carlo 
method. Observe that both expressions (|36ap and <|36b|) 
are written as averages over probability distributions 
P[h(m)] and Q[u(m)] and vectors J. The Monte Carlo 
algorithm that we describe below represents distributions 
P[h{m)\ and Q[u(m)) as finite samples {hi(m)}fL 1 and 
{ui(m)}fL 1 . (Implementation details of storing functions 
h(m) and u(m) in memory are not discussed here; we as- 
sume that it can be done without any loss in precision). 
A single iteration step can be implemented as follows: 
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1. Compute a sample {•Uj(m)}. For each i S 1,...,N 

(a) Choose h2(m), . . . , Kr (m) from the set 
{hi(m)} uniformly at random. 

(b) Choose a disorder vector J at random. 

(c) Evaluate u{m) = uj (m; [{^(m)}^l 2 ] ) using 
Eq. (341). 

2. Compute an updated sample {/i-(m)}. For each 
ie 1,...,N: 

(a) Choose a random integer d from the Poisson 
distribution with parameter Kj. 

(b) If k = 0, let /i'(m) = uo(m), otherwise 

(c) Choose ui(m), . . . ,Ud(m) form the set 
{iti(m)} uniformly at random and 

(d) Evaluate h'(m) using 



h'(m) — uo(m) + Uk(r 



(39) 



fc=i 



The convergence criterion for the algorithm is that step- 
to-step fluctuations are entirely due to the finiteness of 
N, i.e. that both the old {ft.,*(m)} and the updated 
{h'^m)} histograms sample the same probability distri- 
bution. This can be verified by the Kolmogorov-Smirnov 
test for instance [5l| . 

Self-consistency equations and their Monte Carlo im- 
plementation are related to Thouless- Anderson-Palmer 
(TAP) equations discussed in Appendix iBl 



C. Quantum limit (T = 0, T > 0). 

A number of simplifications are possible in this limit. 
Since all integrals over magnetizations have the form 
J dme"^" 1 ', in the limit (3 — > oo they are dominated by 
the minimum value of /(m). The replica free energy func- 
tional TlP[h{m)]\ = jV{P[h{m)]j - S{P[h{m)]j retains 
the form given by Eqs. (|25|) and (|28| . but expressions 
(|26)l and ((29j) for Uj[{hi{m)}] and £{h(m)] simplify to, 
respectively, 



Uj [{ht(m)}} = 



mm 

{m t } 



K 



e=i 



K 



Ej({mi}) + V hdmi) -Vmin[^(m)] (40) 



and 

C[h(m)] = min[/i(TO)], 

m 

while Eq. (|30|l assumes the asymptotic form 



Mo(m) = — r-\/l 



?7Z 



(41) 



(42) 



Self-consistency equations retain the form of Eqs. (|36|l . 
but the expression for uj(m; {{he(m)}]) reduces to the 
following: 



u(m) 



mm 

m,2,...,m K 



K 



Ej(m, m 2 , . . . , m K ) + 2J hi(me) 

K 

Emm\hi(m)]. (43) 



The physical meaning of the effective fields h{m) is par- 
ticularly evident in the limit T — 0. The free energy cor- 
responds to the minimum of the effective Hamiltonian of 
Eq. {3D): 



Wr=o({mi};{J}) 



{ii—in) 



(44) 



In the limit r — > oo the free energy is dominated by the 
second term: F = ~T, which corresponds to a state with 
all spins completely polarized along the x direction. In 
the limit V — > the free energy is expected to be F w 
in the satisfiable phase and F > in the unsatisfiable 
phase. 

For each spin, hi(m) is, up to a constant, the increase 
in energy if the magnetization of spin i is set to m (mag- 
netizations of other spins are allowed to adjust). 

It is possible to set up a deceptively simple system of 
equations for magnetizations {m*} corresponding to the 
minimum of (@4|). Solving dH.T=o/dm,i = we observe 
that m* may be represented in terms of scalar effective 
fields h* via 



h* 



(45) 



while each h* is a sum of contributions Uk from each 
hyperedge incident to vertex i. E.g. for if-SAT 



n 

i=i 



1 + Jem 



(46) 



The description of the problem in terms of order pa- 
rameter P{h*) — the histogram of fields h* — is effective 
for large values of T where (|44|) has only one local min- 
imum. However in the limit of small T the number of 
local minima becomes exponential in N, which is the es- 
sential reason for the introduction of the functional order 
parameter. 



IV. SMALL TRANSVERSE FIELD REGIME AT 
ZERO TEMPERATURE. 

For small values of the transverse field, the free energy 
functional can be expanded in powers of T around T = 
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corresponding to the classical limit. The limit T = has 
been considered in Appendix IA 41 Taking the limit T = 
afterwards, expression J5]) for the classical free energy is 
recovered. We expect that the physically relevant value 
of the free energy (unlike that of the order parameter) 
cannot be affected by the order in which the limits are 
taken. It is instructive to verify that the same result is 
obtained when the limit T = is taken first, followed by 
r = 0. Even though the effective field functions h(m) will 
be finite everywhere in the interval [— 1;+1], the value 
of the free energy will by determined by the values h±\ 



attained on both ends of the interval. 

As a first step, we demonstrate that the function 
uj(m; [{hi(m)}}) is always convex. This convexity prop- 
erty is valid for arbitrary values of T. We evaluate u(m) 
for some linear combination of magnetizations m and 
mi. Writing m%, . . . , m* K to denote the values of magne- 
tization that minimize the first term on right hand side 
of Eq. l(43|) and using the property that Ej(mt, . ■ ■ , ma) 
is a multilinear function of magnetizations, we write the 
lengthy inequality proving the convexity of u(m); 



u(amo + (1 — a)mi) — min 

m 2 ,...,mjc 



A' 



A 



1=2 



Ej(am + (1 - a)rai,TO 2 , . . . , m K ) + h e (mi) - }mm[he(m)] 

i i r ' m. 

1=2 

A A 

aEjOmo, rn%, ■ ■ ■ , rn* K ) + a ) he(m?) + a> mm[he(rn)] 



1=2 



A 



A 



+ (1 - a)Ej(mi,m* 2 , . . . , m* K ) + (1 — a) hAm*p) + (1 — a) mm[hf.(rri)] 

* — ' * — ' m 

1=2 1=2 

^ au(mo) + (1 — a)u(mi). 



(47) 



Usin g the convexity of u(m), it can be established from 
Eq. (|36bip that in the limit r = the effective field 
functions h{m) are also convex due to the vanishing of 
Mo (m). The convexity of h(m) and the multilinearity of 
Ej({ m i}), together, ensure that expressions of the form 
Ej(mi, . . . , tuk) + hi(mi) achieve their minimum values 
for m,£ — ±1. Similarly, minima of effective fields h(m) 
can be replaced by min(/z_i, h + i) due to convexity of 
h{m). It follows that the value of the free energy will be 
unchanged if minima over the interval m G [—1; +1] are 
replaced with minima over the discrete setm 6 {— 1; +1}. 
Hence, the free energy of the quantum model in the limit 
T = must equal that of the classical model. 

One corollary to this is that in the limit T = the 
functions u(m) are piecewise linear. Indeed, du/dm = 
{d / dm)Ej{m, rri^, ■ ■ ■ , rn* K ) may depend on m only in- 
directly via {nig}f =2 . Since m\ G {— 1;+1} the slope of 
u(m) cannot change continuously; instead it assumes one 
of finitely many values depending on the value of m. 

So far we have kept the derivation as general as possi- 
ble. In the following we restrict our attention to random 
if-SAT proper described by the cost function l(27|l . In 
the limit T = functions u{m) (sketched in Fig. ^ may 
be parametrized by a single parameter u as follows: 

u(m) = min (2, 2\u\, 1 — (sgnM)m). (48) 

Using the same letter for the function it (to) and the 
parameter u should not lead to confusion. We will al- 
ways include the magnetization argument to refer to the 
function u(m). The value of u(m) for a particular mag- 



(a) |ii|<l (b) |u|=l 




FIG. 2: Form of u{m) in the classical limit (r = 0). Two 
cases are depicted: (a) |u| < 1, u > and (b) \u\ = 1, u > 
0. Analogous figures for u < may be obtained by mirror 
reflection m — > — m. 

netization (e.g. to = or to = ±1) will be indicated 
using subscripts: i.e. uq, u±i. 

It can be seen from Eq. lf48|) that u — \{u-\ — u + i). 

Although h(m) = uo(m) + J2k=i u k(m) does not admit 
a simple parametrization, we can still define scalar h — 
\{h-\ — h +1 ). This choice ensures that h = J2k=i u k- 
As expected, u(m) defined by Eq. lj43j) assumes the form 
of Eq. I|48p and depends on {hi(m)} only via {hi}: 

u = min(l, ( J2h 2 ) + , . • • , ( Ja^a) + ) , (49) 

with (x)+ used to denote max(a:,0). 

This describes two different regimes. Function 
u{m) has the form depicted in Fig. [2] (left) whenever 
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va\o.i{(Jihi) + \ < 1, and that shown on the right if 
mini{J e (he) + } ^ 1. 

The order parameter P[h(m)] may be obtained by 
iterating Eqs. (|36| starting from, e.g. P^[h{m)\ — 
S[h(m) — Uo(m)] corresponding to the non-interacting 
model [66]. The effects of small, but finite values of T 
can be illustrated by performing a single iteration. Sub- 
stituting ft,2(m) = ■ • ■ = hK{m) = uo(m) into (|43|) gives 



(a) \u\<T 



uj(m) = min 

m,2,...,m K 



1 + Jim 



K 1 

n 1 



Jimt, 



- {K - l)lVl -to 2 



{K-1)T. (50) 



This expression is neither zero (as in the classical case), 
nor even piecewise linear. It should be declared in ad- 
vance that we do not need the precise analytical expres- 
sion for uj (to) as the free energy will not depend on such 
details. It is easily seen that uj(m) is monotonic in to, 
and that it is zero at m — —J\. In addition one can 
demonstrate that 



uj{m) — r — o(r) when 1 + J\m 3> T. 



(51) 



When K ^ 3 this approximate identity is strengthened 
to uj(m) = r for 1 + Jito ^ CT (for some constant 
C). This form of uj(m) is sketched in Figure [3] (left) 
for Ji < (in this particular case u(0) = T). Since we 
are not concerned with the precise form of u(m) it is 
still permissible to describe it using a single parameter 
u = o(ii_i — m+i) (which would equal — JiT/2 in the 
present case). Expression (|48|) would apply everywhere 
on [—1; +1] except for the vicinity of m = sgnw, where 
1 — (sgnu)m = 0(T). Note that if either 1 — (sgnu)TO <C 
T, or 1 — (sgn?i)m >• T, expression (|48| remains valid up 
to o(r). 

By considering additional iterations of Eqs. (|36l) it is 
possible to classify all possible forms of it (to) that can be 
encountered. In addition to the piecewise linear forms of 
Fig. [21 it may have one of the forms depicted in Fig. [3j 
The latter form may occur only if \u\ ^ T/2 (Fig. [21 left) 
or 1 - T/2 < M < 1 (Fig. [21 right). 

Observe that Eq. lj48]) is approximately valid for all to, 
with the possible exception of 1 — \m\ ^> T. Recognizing 
that u(m) is monotonic and that |du/dm| ^ 1, we can 
restate the condition in an equivalent form. We require 
that du/dm approximately [up to o(T)) equal either or 
±1 for 1 — \m\ >> T. For values of to such that 1 — |m| = 
O(l) the derivative du/dm equals either or ±1 with a 
correction of at most 0(T 2 ). 

To investigate the qualitative form of effective fields 
h(m) write Eq. (|39j) substituting the value of uo(to): 



h(m) = —Ty/ 1 — m 2 + Uf.(m). 



(52) 



fc=i 



Function h(m) is a sum of concave and convex functions. 
One of possible forms of h(m) is sketched in Fig. [H All 



4\. 




0(T) ' 

FIG. 3: Possible form of u(m) for finite, but small F. Figures 
depict two possibilities corresponding to u > (u < corre- 
sponds to mirror images rn — > — m): (a) |u| ^ T/2 and (b) 
\u\ > 1 — T/2, For 1 ± m = 0(T), the functions u(m) are not 
piecewise linear. Together with Fig. [2] this encompasses all 
possible forms of u(m) in the limit of small T. 



features that are o(T) have been suppressed. In particu- 
lar, Fig. 2] fails to reflect that the locations of the local 
minima at to = ±1 are shifted by 0(r 2 ) In general, local 




FIG. 4: Typical form of the function h(m), parameterized by 
h and h. In general h — \(h-\ — h+i). The distance from 
the middle minimum at m ~ to the centerpoint of the line 
joining h(— 1) and h(+l) is T — h. All features 0(T 2 ) have 
been suppressed (see discussion in text). 

minima of h(m) away from the endpoints of the inter- 
val [—1; +1] must satisfy dh/dm = 0. Since duk/dm are 
approximately integers for 1 — |m| ^> T, such local mini- 
mum can exist only if J2k dv-k / dm\ m =o « and can only 
be located at m* [up to O(r)]. Neglecting contribu- 
tions of 0(r 2 ) and higher, the free energy is determined 
by values of u(m) and h(m) at to = or to = ±1. 

We will parameterize each of u(m) and h{m) by scalars 
u, u and h, h respectively. We define 



U-i - u+i 



U-i - 2u + u+i 



(53a) 
(53b) 



And h, h parameterizing h{m) of Eq. lj52j) are chosen as 
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follows: 



k—l 
d 

h = y^ufc. 

fe=i 



(54a) 
(54b) 



Note that U and u are not independent variables, but 
are related by 



u = min(|u|, 1 — \u\). 



(55) 



Combining Eqs. (|48j) . (|52|) and l|53p . we obtain for values 
of ft(m) at m = ±1 and m = 0: 



ft±l = ^ \ u k\ ± ft, 



fe=l 



tit- 1 + ft — r. 



(56a) 



(56b) 



fc=i 



The expression ^ fe=1 represents a constant shift 
which must cancel out in the expression for the free en- 
ergy. This cancellation allows one to parameterize ft(m) 
by ft, ft alone. 

It is straightforward to rewrite the self-consistency 
equations lj36|) in terms of the reduced distributions 
P(ft, ft) and Q(u,u). However, it is more instructive 
to derive self-consistency equations from the stationar- 
ity condition for the free energy functional T [P(h, ft)] 
that can be derived by using our ansatz for ft(m). 

As before, we separate the free energy functional into 
two parts corresponding to the quasipotential and the 
quasientropy: T [P(ft, ft)] = 7V [P(ft, ft)] - S [P(ft, ft)] . 
We write down without proof the expression for the 
quasientropy S [P(ft, ft)] : 



S = \ dftdft£(ft,ft) 



dwdw 



\u)h-\-\ujh 



E(u,u), 



where £(ft, ft) and E(ui,oj) are given by, respectively, 

£(ft,ft) =max(|ft|,r-ft), (57a) 



E(u, uj) = P(u, Q)(l- In P(ui, u)j\ 



(57b) 



with P(u,Q) = J dftdft e iu3h+iah P(h, ft) used to denote 
the Fourier transform of P(ft, ft). The derivation of this 
expression is straightforward and relies on the ability to 
replace all minima over magnetizations in the interval 
[— 1; +1] by those over the discrete set m 6 0,±1. 

The derivation of the quasipotential is slightly more 
intricate. The minimum of 



A' 



E' J (m 1 ,...,m K )=2l[ 



1 + J e m e 



1=1 



A 

E 



h t (m t ) (58) 



may occur only for mi, . . . , mjf = 0, ±1. It is unneces- 
sary to consider all 3 K possibilities, however. Let m\ de- 
note the location of a global minimum of hi(m). The lo- 
cation of the global minimum of (|58l) is such that mi = 
or mi = —Jn for some £, while all other magnetiza- 



tions are m^ 



It is never advantageous to have 



more than one magnetization different from m* t as long 
as T < 1. 

Therefore, E'(mi, . . . , mjf) may be written as a mini- 
mum over just K distinct possibilities. After some alge- 
bra we obtain the following expression for the quasipo- 
tential V [P(ft, ft)] : 



1, h e ) x (Uj(h,h))^ 



K 

V = / {dMft^} Y[P(h 
J 1=1 



with Uj(h, ft) given by 



Uj(h, ft) = 2 min { V (Jih e ,r - hi)} , (59a) 



and the definition of ry(ft, e) is 



7?(ft, e) = min (1, (ft) + ) + \{e~ \h\)+ - \{e - \h - 1|)+ 

(59b) 

(the auxiliary function rj(h, e) is sketched in Fig. [5] for 
illustrative purposes). Note that for e ^ 0, Eq. <|59b|) 
reduces to r)(h, e) = min (1, (ft)+). 



v(h,e) 




+ £ 



1 l + e 



FIG. 5: The form of the function r\ = rj(h,e) defined in 
Eq. ([59b). 



It is immediately seen that in the limit T = 0, Eqs. f57j) 
and l|59p . rewritten in terms of P(ft) = J dftP(ft, ft), coin- 
cide with classical T — expressions for the quasientropy 
and the quasipotential respectively. 

The stationarity condition is 6(jV — S)/8P(h,h) = 
const. It should come as no surprise that the following 



13 



identity holds: 

SV _ 
SP(h,h) ~ 

K (^J dudu Q(u, u)C(h + u,h + u) — C(h, h)^j + const, 

(60) 

where Q(u,u) is effectively a distribution of just one pa- 
rameter u: 

Q(u, u) = Q(u)5(u — min(|u|, 1 — M)), (61a) 

K 

Q(u) = / {dhdhe} \\P(hM x u T ;j ({heMf-z) , 
•* i=i 

(61b) 

with u T .j(h 2 , h 2 ;...; h K , h K ) given by 

ur;j({ht, hi}) = e ^ n K {v {Jtht, T — hi)}. (61c) 

Solving the stationarity condition reveals the following 
relationship between P(h,h) and Q(u,u): 



P(h,h)= J dhdhe iu;h+ioh expKjQ(uj,u!), (62) 

where Q(u),Q) is the Fourier transform of Q(u,u). An 
alternative form of (l62l) is 



/d 
{dukduk} Y[Q(u k ,u k ) 

a fc=l 

xs(h-J2 u k)s(h-J2^k). (S3) 



fc=l 



fc=l 



The order parameter can be found by solving Eqs. fBTj) 
and ([62| self-consistently. It is straightforward to write 
down r <C 1 TAP equations for a particular disorder 
realization by "reverse-engineering" these relations, inter- 
preting P(h, h) and Q(u, u) as the histograms of effective 
fields associated with vertices and hyperedges of the ran- 
dom hypergraph. 



V. CLASSICAL ZERO-TEMPERATURE 
SOLUTION REVISITED 

A. Scale-invariant replica-symmetric solution 

The analysis of Ref. [2l| presents a classical picture of 
the first-order phase transition: as a competition between 
two locally stable solutions: the trivial P(h) = 6(h) and 
the non-trivial P{h), However, a Monte Carlo study re- 
veals that the non-trivial solution is not stable for any 
7 < 7 C w 4.60 under iterations of the self-consistency 
equations for P{h) and Q(h). This casts doubt on the 



picture of competition between two local maxima (the 
free energy must be maximized) or the prediction of the 
dynamic transition at jd ~ 4.43. We claim that the phase 
transition at 7 C w 4.60 is in fact continuous. While the 
value of 7 C has been determined correctly, the discontinu- 
ity of the order parameter is an artifact of the discretiza- 
tion used in the numerical procedure (values of (Ah) -1 
up to 30 have been used in [21(). 

In this section we consider the model described by the 
free energy functional (JU, but with the expression J9]) 
modified to 



U 



(O) 



(hi,...,h K ) = 2 min {( J e h e ) + } . 
i=i,...,K 



(63) 



We will refer to this modified version as Model 0. The 
original version will be called Model A. 

The distinguishing feature of Model O is the absence 
of any explicit scale. The free energy functional becomes 
covariant with respect to scaling transformation (rescal- 
ing of effective fields by a factor of A) : 



T [ - 0) [XP{h/\)]=\T {0) [P{h)]. 



(64) 



The immediate consequence is that the maximum value 
of [P(h)] can be either or +oo, depending on the 
value of 7. We can still formally write self-consistency 
equations satisfied by P{h) and Q(u): 



K 

Q(u) = J {dh e } [J P(h e ) x *( u - uj (h 2 , h K )), 



(65a) 



P(h) - E MK-f) f {du k } [J Q(u k ) x 5 (h - «*) , 

d J k=l ^ k ' 

(65b) 

however Uj(h,2, • ■ ■ , hx) becomes linear in I12, ■ ■ ■ , ha'- 
u { °\h 2 , ■ ■ -,h K ) = -Ji f _m in K {{Jihi) + } ■ (66) 

Under successive iterations of self-consistency equations 
([66"1 the distribution quickly converges to a universal 
form, with any subsequent iterations merely rescaling ef- 
fective fields by a factor of A that depends on the value 
of 7: 



ptr+V(h) = XP (r) (h/X). 



(67) 



It is convenient to introduce the simplified order param- 
eter: the width A of distribution P(h). One possible 
choice for the definition of A is 



dhP(h)\h\. 



Successive iterations rescale the value of A by a factor of 
A so that it flows towards one of two fixed points: A* = 
or A* = +00, depending on the value of A. Since A(7) is 
expected to be monotonically increasing, we expect that 
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A* = whenever 7 < j c and A* 
value of 7 C being a solution to 

A(7c) = 1. 



-00 for 7 > 7 C ; the 



(69) 



It is easy to verify that the condition A = 1 is equiva- 
lent to T[P*(h)] — 0, where the free energy is evaluated 
for the limiting distribution P*(h) rescaled so that it has 
finite non-zero weight. This universal form of P(h) con- 
tains a delta-function peak at ft. = as well as a contin- 
uous part. 



(70) 



P(h) = (l-q)5(h)+qp(h), 



where p(h) has been normalized to unity. From Eq. i|66j) 
we can deduce the general form of Q(u), which we sepa- 
rate into a delta function peak at u = and a continuous 
part x(w): 



Q{u) 



K-l 



8{u) 



K-l 



X{u). (71) 



Because an iterated convolution of x( u ) ls a l so continuous 
we can write self-consistency equation for the singular 
part of P{h) and Q(u) only: 



(72) 



A similar equation appears in the analysis of the leaf- 
removal algorithm for random if-XOR-SAT (H- The 
correspondence becomes exact with the replacement 
-f/2 K ~ 1 — > 7. Equation lj72j) admits two solutions: a triv- 
ial solution q = as well as a non-trivial solution q > 
which appears discontinuously above some threshold. For 
K = 3 the critical value of 7 is -y q « 4 x 0.818 « 3.272. 
This threshold is irrelevant for our problem, because the 
corresponding A(7 g ) < 1 and P(h) = 5(h) maximizes the 
free energy. 

Equation lj72j) has the following interpretation. We 
identify q with the fraction of almost frozen variables: 
variables that take the same value for all configurations 
with the lowest energy except for an exponentially small 
fraction. We randomly choose a spin variable and the 
corresponding vertex in associated hypergraph (call it a 
cavity vertex). The degree of this vertex (the number 
of hyperedges incident to it) is Poisson-distributed with 
mean Kj. Each hyperedge connects the cavity vertex 
to K — 1 neighbors (see Fig. [6]). Each of these vertices 
corresponds to another almost frozen spin with probabil- 
ity q; since we expect that spins are equally likely to be 
frozen to +1 and —1, with probability (q/2) K ~ 1 the ef- 
fect of the corresponding constraint is to force the cavity 
spin to have a value of —J\ in all but an exponentially 
small fraction of configurations with the lowest energy. 
The cavity spin will be almost frozen if this happens for 
at least one hyperedge. Since the cavity spin has been 
chosen randomly this probability equals q, which leads to 
the self-consistency condition (|72|) . Observe that our in- 
sistence on variable being almost frozen rather than com- 
pletely frozen is crucial. It may happen that two or more 




FIG. 6: Hypergraph with K = 3. Cavity vertex (white circle) 
has d = 3 (as shown in this picture) hyperedges incident to 
it. Each hyperedge connects the cavity vertex to K — 1 = 2 
neighboring vertices (black circles). In the absence of a cav- 
ity vertex (and incident hyperedges), each of the neighboring 
vertices would be almost frozen with probability q/2. The 
cavity vertex is almost frozen with probability q. The self- 
consistency condition on q that takes into account the Poisson 
distribution of degrees d is given by Eq. (1 72 1) . 



constraints satisfy the condition described above, and ex- 
actly half of constraints have J\ = +1 while the other half 
have J\ = — 1. The net effect is that the cavity variable 
remains unfrozen as a result. When we speak of almost 
frozen variables, we assume that in the thermodynamic 
limit it is unlikely that an exact cancellation takes place, 
i.e. that the number of configurations with sq = +1 and 
so = — 1 is roughly the same (neither is exponentially 
smaller than the other). This simplified picture should 
work if the number of completely frozen variables is much 
smaller than the number of almost frozen variables. 

Once the nontrivial solution q > to (|72|) is found 
(there will be two solutions q > 0, but only the larger one 
is stable) , we can write the self-consistency equations for 
the continuous parts p(h) and x( u ) : 



X (U) = (K-I)p(u)[ J <\hn(h) 

d 

T 



K-2 



(73a) 



d=l 



fe=i 



(73b) 



(the "renormalized" connectivity is 7 = K^q K 1 ). These 
equations can be solved iteratively. The value of A is 
estimated after each iteration by rescaling p(h) so that 
J dhp(h)\h\ = 1. The critical value 7 C is found from 
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A(7c) = 1- We obtain j c « 3.1650 which translates into 
7 C f» 4.6002. This agrees with the threshold value re- 
ported in Ref. [ll|. In Fig. [7] we plot p*(h) at the critical 
value of 7. 



0.5 



0.4 



0.3 



p(h) 



0.2 



0.1 




h 



FIG. 7: Scale-free solution p*(h) of equations l|73p . for 7 = 7 C . 
Inset shows In to illustrate the approximately exponential 
decrease of p(h) as h — * 00. 



B. Continuous phase transition in classical T = 
K-SAT 



Armed with the solution to Model obtained in the 
previous section, we can make qualitative predictions 
about the solution to classical T = A'-SAT (Model A). 
The self-consistency equations for Model and Model A 
differ only in the expression for Uj{h%, . . . , hn): 



u 



(A) 



(h 2 



,h K ) = -J x min(l, 



,(0) 



(l>,2 



,h K )\ 



(74) 

In the limit A -C 1 both expressions are nearly equal, 
A defines the scale of /12, • • • , fiK- As a result, 



since 



we expect that whenever 7 < 7c , successive iterations 
of self-consistency equations for Model A flow toward 



0. The fact that 



,(^)| 



< u 



(O)l 



erate the process. Conversely, when 7 



that 



,( A )\ 



will only accel- 
> 7c , the fact 



1 j I < 1 prevents the divergence of A, which will 

stabilize at A = 0(1). The value of 7 = 7c°' ) is the 
boundary between the satisfiable (A = 0) and the unsat- 
isfiable (A > 0) phases in Model A. 

The dependence of the order parameter A on connec- 
tivity 7 in the vicinity of phase transition I7— 7 C | <C 1 can 
be estimated variationally. We expect that A increases 
continuously from the value of at 7 = 7 C . Right above 
the transition, in the limit A <C 1, Model O and Model A 
are essentially equivalent. For the variational ansatz for 
P(h), we choose the scale- invariant solution of Model O, 
P A (h), corresponding to 7 = 7 C . The width A of the 
distribution appears explicitly and is the adjustable pa- 
rameter. Exactly at 7 = 7 C , the free energy of Model O 



is degenerate (_F(°) = 0). For Model A, this degeneracy 
is lifted, and we can obtain the width of distribution A 
by optimizing 



K^(A) 



dh 



dh'P A (ti) 



K 



V ( °\P A (h)} 



dh 



dtiPKh'] 



K 



(75) 



The quasientropy is independent of the choice of the 
model: S${A) = S (0) [P A (h)}. Observe that 



s^[P* A (h)] = icV^[PUh)}- 



(76) 



The asymptotic form of P A (h) is related to that of p*{h) 
(see Fig. [TJ inset): 



p*(h) ae-" (h)|h| , 



(77) 



where p(h) is a function of very slow growth. In par- 
ticular, it grows slower than iterated logarithm of \h\. 
Therefore, in the limit A « 1, the correction term in 
(HH scales as e '^^ A ^ A . 

The variational free energy may be written as follows: 



A var (A) « a( 7 - 7c) A - e -^ A ^ A . 
Solving di^ /dA = with respect to A yields 

1 



7- 



7c « ^2 e 



- M (l/A)/A 



(78) 



(79) 



With some abuse of notation (we write x ~ y to mean 
that x is asymptotically proportional to y with the coeffi- 
cient of proportionality being an extremely slow-varying 
function of y) , the dependence of the order parameter A 
on connectivity 7 > j c may be written as follows: 



1 



|ln(7-7c)| 



(80) 



Given the extremely singular character of this function, it 
is not surprising that the transition looks like a first order 
transition in numerical simulations. Critical exponents 
a = 1, (3 = are precisely those expected for the first 
order transition (the scaling exponent associated with the 
logarithm is zero). In the vicinity of the phase transition, 
just above it, the behavior of the free energy is 



7 -7c 



|ln(7-7c)| 



(81) 



Let us briefly describe the mechanism of this contin- 
uous phase transition. In the classical case, the free en- 
ergy is the maximum of the free energy functional [over 
symmetric distributions P(h)) as demonstrated in Ap- 
pendix IA4I 



F( 7 )=max[7V[P(/ i )]-5[P 



(82) 
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It is convenient to regard the reciprocal of the connectiv- 
ity I/7 as a Lagrange multiplier for S[P(h)] . We consider 
the quasipotential V as a formal function of the quasien- 
tropy S: 

V(S) = m^[V[P(h)]\S[P(h)] = Sj. (83) 

P(h) 

In Fig. [8] we sketch the function V(S), which should be 
convex. 



A 

tan a = 7c 




FIG. 8: Quasipotential as a funtion of quasientropy in Model 
A. Tangentials to the curve have local slope 7 _1 . There are 
no tangentials with slopes corresponding to 7 < 7c , hence 
S — V = in this region. 

Via newly defined V(S), the free energy F( 7 ) can be 
expressed as follows: 

F( 7 )=max[ 7 V(S)-S]. (84) 

The value of S that maximizes the right hand side of |84|) 
is a solution to 

7- 1 = dV(S)/dS, (85) 

for 7 7 c . When 7 < 7c , Eq. (|85|) has no solutions and 
the r.h.s. of f84|) is maximized by S — V = 0. 

Iterating self-consistency equations <[65|l for vari- 
ous values of connectivity 7 becomes extremely time- 
consuming in the vicinity of the phase transition. We 
have verified the general trend that iterations converge 
to A = for 7 < 7 C and to a finite value of A for 7 > 7c . 
To find the equation for the order parameter A( 7 ) nu- 
merically, we took a different route. Instead of fixing 
the value of 7 and iterating equations until convergence, 
we fix the width A of probability distribution P(h) and 
choose the value of 7 at each iteration step so that P(h) 
has the desired width A. In fact, by making an appro- 
priate choice for the somewhat arbitrary definition of A, 
the corresponding value of 7 may be obtained at every 
step with just a single arithmetical operation. One such 
choice — A = J ' dh P(h)h? — exploits the identity 

J dhP(h)h 2 = K~/J duQ(u)u 2 . (86) 

The above-described approach results in tremendous 
speed-up. While the number of iterations required for 



convergence for fixed 7 increases to infinity as 7 ap- 
proaches 7c, for fixed A a complete convergence (to ma- 
chine precision limit) is achieved within 20 iterations. 

Additionally, we utilize a Quasi Monte Carlo (QMC) 
method [H3| that we have formulated specifically for 
the problems involving probability distributions. Un- 
der ideal conditions, the expected error for Quasi Monte 
Carlo (QMC) is 0(logN/N) compared to 0(l/y/N) for 
standard Monte Carlo (MC). In practice, discontinuities 
and singularities worsen the error estimate, however the 
asymptotic behavior of QMC is always better. For values 
of N as small as N = 8192 the error in the value 7 ob- 
tained using standard Monte Carlo was 1%, while Quasi 
Monte Carlo produced an error of 0.05%. The advantage 
of QMC over MC only increases with increasing N. We 
have used values of TV up to N = 2 25 = 33 5 5 4 4 32. An 
overview of our method is presented in Appendix O 
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FIG. 9: Connectivity 7 vs. width of the distribution of ef- 
fective fields A in Model A. We predict that 7 > 7 C for any 
A > 0, so A(7) has no discontinuities. The inset replots the 
data to illustrate the asymptotic relation J79}. The y axis 
corresponds to /( 7 , A) = [3.5 — In A 2 (7 — 7c )] _ , which is 
asymptotically linear in A as A — > 0. Dotted line is the re- 
sult of extrapolation to small values of (7 — 7c) where the 
numerical error is too large. 

In Fig. [9j we plot the function 7(A) obtained numer- 
ically. To establish that the phase transition is continu- 
ous, we must convince ourselves that 7(A) is a strictly 
increasing function of A, i.e. that 7(A) > 7(0) for ar- 
bitrarily small A. The inset shows the roughly linear 
dependence of /( 7 , A) = [C - In A 2 ( 7 - 7c)]" 1 on A (cf. 
Eq. [79j) that, when extrapolated, predicts the vanishing 
of A as 7 — > 7c. 

In Fig. [lOl we plot the dependence of the free energy 
F on connectivity 7 in the vicinity of j c . Contrary to 
visual perception the slope of F(j) at 7 = 7 C + is zero 
from Eq. (|8l)l . The "apparent" slope AF/A7 decreases 
as a function of A as can be seen by comparing the main 
figure with the inset in Fig. [TU1 
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FIG. 10: Free energy F vs. connectivity 7. Contrary to visual 
perception, F(fy) is not linear for 7 > j c . The inset zooms 
in on the transition, keeping the aspect ratio the same. The 
apparent slope of ^(7) in the inset is smaller. This apparent 
slope will tend to zero as progressively higher zoom ratios are 
used. 



VI. NUMERICAL RESULTS 

A. Classical regime (r = 0, T > 0) 

In addition to Model O and Model A described in 
Sec.[Vl we introduce two new models: Classical Model B 
and Classical Model AB. Classical Model AB is precisely 
the finite temperature classical random if -SAT. The dis- 
tinguishing feature of Model B is the absence of explicit 
temperature. It is defined using Eqs. (|25|) . (|28| . but with 
the following choice for Uj{{h{\) and C(h) respectively: 



^ (i + e -2Jihi )... (1 + c -2J K h K ) J ' 

(87a) 

C {B) =ln(2coshft). (87b) 

All four models (0,A,B and AB) can be described by a 
single form of the free energy functional that depends ex- 
plicitly on two parameters: the temperature T = 1/(3 and 
the energy scale parameter A. This common model can 
be defined using the following expression for Uj{{hf}) 
and C(h): 

( 1 - e" 2A / T \ 
Wr A = -Tln 1 — , (88a) 

v nf =1 (i+e- 2 w)y 

£ t ,a = T In ^2 cosh ~ \ . (88b) 

We summarize the values of T and A for the four models 
we have introduced in Table |H 



TABLE I: Four different models defined by the values taken 
by parameters T, A. Statistical properties of Models O, A, 
and B should be similar, since T/A = in all three cases. 
Model AB is the classical finite-temperature random Tf-SAT. 



Type of Model 


Temperature 


Scale A 


Model O 


T = 


A = oo 


Model A 


T = 


A = 1 


Model B (classical) 


T = 1 


A = oo 


Model AB (classical) 


T > 


A = 1 


This common model with explicit dependence 


on T,A 


satisfies the following 


scaling relations: 




FT,A[P(h)} 


= TT T =iATP{h/T)l 


(89a) 




= A^ t ,a=i[AP(/i/A)]. 


(89b) 



The implication is that the statistical-mechanical proper- 
ties of this model depend on the ratio of two scales T/A. 
In particular, we expect that Model O, Model A and 
Model B undergo phase transition at the same value of 
the critical connectivity, since T/A = in all three mod- 
els. We have previously established that = 7e°' ) . 

The numerical results for the Classical Model B are 
presented in Figs. [TT] and [321 To obtain the numerical 
solution, we adopted the same strategy as for Model A. 
We computed 7 as well as a number of other quantities for 
each value of A. An interesting feature of Model B is that 
the function 7(A) plotted in Fig.[TT]is non-monotonic and 
cannot be inverted unambiguously for 7 > 7 < ^ s - ) (+oo). 
Although not reflected in the figure; formally there exists 
another solution corresponding to A = +00, with the 
free energy F = +oo(— 00) for 7 > 7a (< l^)- Since a 
branch with the higher free energy must be chosen (see 
Appendix IA 4[) . the branch A* < A < +00 is unstable 

and T( s ) = +00 for 7 > 7c°^. This behavior may be 
understood in terms of the form of function V(S). The 
discontinuous transition occurs because it has non-convex 
form as sketched in Fig. fT3l 

Whereas the free energy of Model A (see Sec. EE) 
corresponds to the internal energy E of Eq. J7]) (or 
the number of violated constraints); the free energy of 
Model B is related to entropy S: 

- NFW = <E) = (InAfs), (90) 

where Afs represents the number of solutions that sat- 
isfy all constraints. The divergence of F^ signals the 
transition to the unsatisfiable phase (Afs = 0). 

In Fig. Q2I we plot the specific entropy (i.e. the negative 
of the free energy) as a function of A. Note that since 
the entropy is finite at 7 = 7 C , the number of solutions, 
just prior to the satisfiability transition, is exponentially 
large. It is the expected behavior: for the associated 
hypergraph, random graph theory [H3| predicts that there 
are O(N) vertices that are either isolated or belong to 
small isolated clusters. These make a finite contribution 



18 



4.62 



4.60 - 



4.58 



4.56 



4.54 



4.52 








- 7c f 

1 . 


A* 

i 



io~ 2 io _1 10° io 1 Aio 2 io 3 io' 



10 



A 



100 



FIG. 11: Connectivity 7 vs. width of distribution of effective 
fields A in Classical Model B. The main figure shows 7(A) in 
the region close to critical. The inset shows 7(A) in a wider 
range. The branch A > A* is unstable, hence A = +00 as 
soon as 7 > 7c- The unstable branch is drawn with the gray 
solid line. 
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FIG. 12: Specific entropy S/iV (logarithm of the number of 
solutions) as a function of A. The main figure shows the 
dependence in the critical region, while a wider range of A is 
used for the inset. The unstable branch (A > A*) is drawn 
in gray. The entropy decreases to E* as 7 approaches j c , 
but jumps to -co (corresponding to zero solutions) for 7 > 
7c. Entropy corresponding to very large values of A in the 
unstable branch could not be determined with good precision. 
Results of extrapolation are indicated using the dotted gray 
line. 



to the entropy but do not affect the overall satisfiability 
of the random instance. 

Based on results for Model B we expect that the non- 
monotonic behavior of 7(A) persists for some sufficiently 
small but finite temperatures. In Fig. [14] we plot the 
functions 7(A) for the classical Model AB for a range of 
temperatures from T = 0.01 to T = 0.5. It is seen that 
far away from 7 = 7 C Model AB interpolates between 
the regimes of Model B with A = 0(T) and Model A 
with A = O(l). For small temperatures 7(A) is non- 




FIG. 13: Illustration of the non-convex behavior of V(S) for 
Classical Model B. Slopes of tangentials to the curve deter- 
mine values of 7 via Eq. lf85j) . The largest attainable values 
of the quasipotential and the quasientropy, V* and S* corre- 
spond to 7 = 7c- 



monotonic, which gives rise to the first-order phase tran- 
sition. 




FIG. 14: The dependence 7(A) in Classical Model AB for a 
range of temperatures. Curves labeled (l)-(6) in the inset 
correspond to temperatures T = 0.01, 0.02, 0.05, 0.1, 0.2, 
and 0.5 respectively. The temperatures are labeled explicitly 
in the main figure. The curves smoothly interpolate between 
the regime of Model B [A = 0(T)\ and that of Model A 
[A = O(l)]. The inset shows that 7(A) is not monotonic for 
sufficiently small temperatures. 

We make two dimensional parametric plots 
(7(A), F(A)) for a range of temperatures T (see 
Fig. [T5|. For T < T* these curves are self- intersecting. 
Stable branches (black solid lines) have a discontinuous 
slope at the point of self-intersection which leads to the 
discontinuity of the order parameter. The dashed green 
line (marked (2) in the figure) is the line of singularities 
between (7c, 0) and (j*,F*) terminating at the critical 
point. In the space of variables (7, T, F) stable and 
unstable branches of F form a dovetail singularity. It 
should be recalled that for T — the derivative dF/d7 
has no discontinuity, although it is difficult to see from 
the figure. 

Finally, in Fig. [T6l we show the numerical phase dia- 



19 




FIG. 15: Two-dimensional parametric plots (7(A), F(A)) for 
the Classical Model AB (finite temperature K-SAT). Differ- 
ent lines correspond to different temperatures. Black solid 
lines correspond to the stable branches of the free energy; 
gray lines correspond to unstable solutions. Switching be- 
tween stable branches occurs along the green dashed line (2). 
Along this line, the first derivative dF/dj of the free energy 
has a discontinuity. Red dashed lines (1) and (3) are the 
spinodals dV/dS = 0. Points A and B along T = line 
correspond, respectively, to the critical threshold in Model A 
and the metastable solution 7 w 4.6184 of Model B. Lines 
(1),(2), and (3) meet at a critical point C, corresponding to 
T* ~ 0.05864. Note that the line BC is nearly vertical, with 
corrections that are ~ exp(— 1/T). 



B. Quantum regime (F > 0, T = 0) 

We introduce Quantum Model B and Quantum 
Model AB as follows. We will keep the definition lj57|) 
for the quasientropy, but in the definition (f59|) for the 
quasipotential the function r)(h, e) will be replaced with 

r] A (h, e) = min (A, (h)+) + l(e - \h\)+ - i(e -\h- A|)+ 

(91) 

so that the free energy functional will contain a character- 
istic scale of the effective fields A explicitly. By choosing 
the values of T and A according to Table [II] we define 
the two quantum models: Model B and Model AB. The 
purely classical models — Model O and Model A — cor- 
respond to the limit T = 0. 

TABLE II: Four different models defined by the values taken 
by parameters F, A. Statistical properties of Models O, A, 
and B should be similar, since T/A — in all three cases. 
Model AB is quantum random A"-SAT; we study the limit 

r < 1. 



Type of Model 


Transverse Field 


Scale A 


Model O 


r = 


A = 00 


Model A 


r = 


A = 1 


Model B (quantum) 


r = 1 


A = 00 


Model AB (quantum) 


r > 


A = 1 



gram in the plane (7, T), The discontinuity of the order 
parameter becomes zero at both ends of the phase bound- 
ary between (p/ c , 0) and {j*,T*), 
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FIG. 16: Numerical phase digram of the Classical Model AB 
(finite temperature A"-SAT). The phase boundary (phase 
transition line) starts from T = 0, critical connectivity 7 C 
and terminates at the critical point with 7* « 4.6185 and 
T* ?s 0.05864. 



Ordinarily, the free energy of the quantum model cor- 
responds to the smallest eigenvalue of the Hamiltonian. 
However, in the limit A = +00, where A defines the en- 
ergy scale of the classical Hamiltonian, the contribution 
from the states with energy E > vanishes. The value 
of the free energy for an instance of Model B may 
be evaluated using the degenerate perturbation theory. 
In the limit A = +00, this free energy is proportional 
to r, which can also be seen from scaling analysis. We 
choose r = 1; the free energy is directly related to 
a property of the space of solutions cr. 

Consider a graph Q having Ms vertices corresponding 
to spin configurations that satisfy all constraints. We 
draw edges between vertices of Q corresponding to con- 
figurations <t,<x' that differ by a single spin-flip. Let A 
denote the adjacency matrix for this graph, i.e. 



A, 



1 if d(tr,tr f ) = 1, 

if d(tr, a') = or d(cr, cr') > 2, 



(92) 



where d(cr,cr') denotes the Hamming distance between 
spin configurations cr and cr' . The free energy of Model B 
will be related to the norm of matrix A: 



nfw = (Mil) - (A raax (^ CTCT ')> 



(93) 



(the spectrum of -4 CT<T ' is symmetric). 

The expression for the free energy may be simplified 
in the limit A = +00. Since in this limit, u = \u\, and 
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u is equally likely to be positive or negative, we may 
express joint probability distribution Q{u,u) in terms of 
the probability distribution of u denoted Q+(u): 

Q(u,u) = ^(Q+(u)S(u - u) + Q+(u)6(u + u)). (94) 

Substituting this into Eq. (|62| . we ontain the following 
factorization of the joint distribution P{h, h): 



. — . fh + h\ (h — h 

P{h,h)=P + (—) P+ (— 



P+(h) = I dhe^exp— ^ 



(<5+(w) 



i 



(95a) 
(95b) 



We can use relations (f94|l and (|95|) to write the free energy 
as a functional of P+(h) alone. The resulting expression 
for F[P + {h)} is 

K 

T^[p+{h)] = 7 / {d m } J] Rfa) x 2 fc min K {^i 

** IS 1 1 " ' ' 



d/iidft,2 C(hi, h 2 ) 



dwidw2 
2tt 



iwihi+ki^Tis 



where the expressions for £(hi,h 2 ) and U(u;i,a;2) are 

h 2 ) = max(|/ii — /i2|>r — h), (96a) 

£(o;i,£J2) = P+(w 1 )P + (w3) (l-lnP + (cj 1 )-ln J P + (w 2 )) 

(96b) 

[-P(w), as usual, denotes the Fourier transform of -P(/i)]- 
The distribution R(rj) that enters on the r.h.s. of 
Eq. f96|) may be related to P+{h) as follows: 



flfa) = J dh 1 dh 2 P+(h 1 )P + (h 2 ) 

x £(77- (max(r/2,/n) 



The numerical results for Quantum Model B are pre- 
sented in Figures [17] and [18l In contrast to Classical 
Model B, 7(A) is a monotonically increasing function of 
A (See Fig. [lTJ) . Its inverse A (7) is a single- valued func- 
tion exhibiting no discontinuities. It diverges as 7 ap- 
proaches 7 C . It is fortunate that there is a single branch, 
as the stability analysis is more complicated in the quan- 
tum case. 

In Fig. [HI we plot ||.A||/iV: the norm of the matrix de- 
scribing the connectivity of solutions. It is seen that this 
quantity does not go to zero as 7 — > 7 C . This can be ex- 
plained by the effect of small clusters in a random hyper- 
graph associated with an instance. This hypergraph is a 
collection of isolated clusters: a giant cluster of size O(N) 
and a large [O(iV)] number of small clusters. Each clus- 
ter may be used to define a subinstance of the problem. 
The space of solution of the large instance is a Cartesian 
product of spaces of solutions of subinstances. It can be 



shown that the norm for the full instance may be 
written as a sum of norms ||*4fc|| for all the subinstances 
corresponding to isolated clusters. The large number of 
small clusters contributes to the finite value of ||«4|| as 
7 -> 7c- 
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FIG. 17: Connectivity 7 vs. width of distribution of effective 
fields A in Quantum Model B. The main figure shows 7(A) 
in the region close to critical. Since 7(A) is monotonic, there 
is only a stable branch. The inset shows 7(A) in a wider 
range. As 7 approaches y c the value of A increases to infinity 
conntinuously. 
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FIG. 18: The norm (largest eigenvalue) of matrix A, defined 
by Eq. l[92|l . as a function of A. The main figure shows the 
dependence in the critical region, while a wider range of A is 
used for the inset. There is only a stable branch. Values of 
||^4|| that could not be computed reliably for very large A are 
estimated using extrapolation, which is indicated by the use 
of black dotted line. 

We should mention that the computed value of 
is not quantitatively correct even in the regime where 
the replica-symmetric solution is stable. This is due to 
our making a static approximation. Although Quantum 
Model B describes the limit V — > 0, the static approxima- 
tion requires a stronger condition j3T — > in order to be 
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exact. We, however, work in the opposite limit j3T — > oo. 



A 




FIG. 19: The dependence 7(A) in Quantum Model AB for 
a range of transverse magnetic fields. Curves labeled (l)-(7) 
in the inset correspond to transverse fields Y = 0.001, 0.002, 
0.005, 0.01, 0.02, 0.05, and 0.1 respectively. The curves in 
the main figure are explicitly labeled with values of Y. The 
curves smoothly interpolate between the regime of Model B 
[A = 0{Y)\ and that of Model A [A = 0{l)\. In contrast to 
the classical case, functions 7(A) are monotonic and free of 
singularities. 

Numerical results for Quantum Model AB are pre- 
sented in Fig. Qjjl We plot 7(A) for transverse field T 
ranging from T = 0.001 to T = 0.1. It can be seen 
that functions 7(A) are always monotonic. In contrast 
to Classical Model AB, the free energy does not exhibit 
non-analytic behavior. The continuous phase transition 
present for T = disappears and is instead replaced by a 
smooth crossover for arbitrarily small T > as depicted 
in Fig. [2Ql The effect of the critical point (7 = -f c , T = 0) 
is that the width of the transition A7 goes to zero to- 
gether with r. We conjecture, by analogy with quantum 
phase transitions in physical systems, that the character- 
istic width of the transition scales as some power of T, 



A7 oc r 



l/z 



(98) 



where the width of the transition has been formally de- 
fined as follows: 



A7 



max 

7 



1 d 7 
AdA 



(99) 



This power law may be verified by plotting points (A7 
and F) on a log-log plot (see Fig. [2lj) , For small T, the 
data seems to converge to power-law scaling with scaling 
exponent z = 1 (the slope corresponding to z = 1 is 
indicated with the gray solid line). However, we have not 
studied this scaling dependence analytically and cannot 
completely rule out the possibility that the dependence 
of the width of the transition on T is more complex and 
cannot be described by a simple power law. 
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FIG. 20: Illustration of crossover transition for quantum PC- 
SAT. The sharp phase transition predicted in classical -K-SAT 
is the critical point at Y — 0. For small but non-zero values of 
F it is replaced by a smooth crossover transition of finite width 
between the underconstrained (7 < j c ) regime described by 
Model B and the overconstrained (7 > 7c) regime described 
by Model A. The width of the critical region decreases as 

r^o. 
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FIG. 21: Transverse field Y vs. the width of the transition 
A7. Log-log scale is used to obtain a power-law fit between 
the width of the transition and the transverse field Y = (A7) z . 
Filled triangles correspond to numerical estimates of A7 for 
the following values of Y: 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 
and 0.1. The slope of the gray solid line corresponds to 2 = 1. 



VII. CONCLUSION 

The main result of this paper is that the thermody- 
namic phase transition between SAT and UNSAT phases 
in classical random i^-Satisfiability problem does not sur- 
vive when quantum effects are incorporated via coupling 
to the external transverse magnetic field. We have stud- 
ied the free energy as a function of connectivity 7 for 
different values of transverse field T. The case T = cor- 
responds to the purely classical limit, and there exists a 
phase transition when 7 is crossing its critical value. We 
have demonstrated that for any small value of T the free 
energy becomes analytic and the sharp phase transition 
at r = is replaced by a smooth crossover transition. 
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This stands in contrast to classical if-Satisfiability model 
at finite temperatures, where we have found a first-order 
phase transition line on the temperature-connectivity 
plane approaching continuously a zero temperature limit. 
However, it is not inconceivable that the seeming differ- 
ence between the classical and quantum case is an arti- 
fact of the replica-symmetric approximation. The RSB 
analysis of dilute antiferromagnetic Potts glass at finite 
temperature indicates that analogous zero-temperature 
static transition becomes a smooth crossover at finite 
temperature [H3|. Whether the inclusion of RSB in T > 
classical if -SAT will similarly lead to the smoothing of 
the static transition is open to investigation. 

We believe the above-mentioned phenomenon is not 
universal among dilute long range spin glasses. We ex- 
pect that in models with Ising or p-spin interactions the 
phase transition at 7 = j c is not affected by trans- 
verse fields below a certain threshold (T < T c ) and that 
the phase boundary has the form of line labeled (1) in 
Fig. [TJ We attribute the difference to the fact that con- 
straints involved in Viana-Bray and dilute p-spin models 
are "stronger" (i.e. a greater number of spin combinations 
are penalized) than those in if -Satisfiability. 

The limitations of our approach are the assumptions of 
the replica symmetry and the failure to include the time- 
dependence of correlation functions. The latter approx- 
imation has been justified on the ground that we work 
in the limit of small transverse fields, as we establish the 
absence of the phase transition line on (7, T) plane. It 
is known that the replica-symmetric approximation can 
capture the existence of the thermodynamic transition 
in the classical if-Satisfiability model while providing an 
overestimated value for the transition point (critical con- 
nectivity) . 

Previously, we attempted to analyze the 0(T) correc- 
tions to the free energy of the r = (classical) if -SAT 
problem along similar lines (HfjJ. In contrast to the 
present analysis, we computed corrections to the integer- 
delta-peaks solution lfl2|) of zero-temperature classical 
if -SAT developed in Ref. plj . Although the assump- 
tion of integer values of effective fields is more natural 
for T = 0, T = and has been used to construct an RSB 
theory of T = K-SAT at the replica-symmetric 
level it does not give a truly stable solution for either 
T > or T > 0. At the same time, the solution derived 
in the present paper is globally stable in the RS sector 
at finite temperatures and transverse fields. 

To obtain a correct location of the phase transition one 
needs to take into account the spontaneous breaking of 
the replica symmetry However, we hope that our 

main result — the smoothing of the phase transition at 
finite values of the transverse field — will be immune 
to the effects of the replica symmetry breaking. On the 
other hand, it is well known that the replica symmetric 
approximation fails to account for the dynamic transition 
and the complex structure of local minima existing in 
the classical if -Satisfiability model at connectivity that 
is smaller than that of the static transition. From the 



perspective of quantum adiabatic algorithm, the likely 
conclusion is that the bottleneck of QAA may be in this 
dynamic transition rather than the static transition. Re- 
cent results on K-SAT show the presence of another tran- 
sition: so-called condensation transition (which coincides 
with the dynamic transition for K = 3) [HIHE]- It is be- 
lieved that the crossover to the exponential complexity 
of classical algoritm happens at the condensation tran- 
sition. It may also be relevant for the performance of 
quantum adiabatic family of algorithms. 

Since the non-analytic behavior of the free energy as- 
sociated with the static transition is the isolated singu- 
larity, it should be irrelevant to the complexity of QAA 
except when 7 = 7 C precisely. It is conceivable that the 
complexity of QAA will be subexponential for 7 7^ j c if 
singularities of the free energy associated with e.g. con- 
densation transition are weaker than that of static (sat- 
isfiability) transition. 

The analysis of a quantum version of if-Satisfiability 
problem lead us to the re-examination of the static phase 
transition in zero-temperature classical limit. In Ref. [2l| 
it has been predicted, using a replica-symmetric analy- 
sis, that this phase transition is of a random first-order 
(discontinuous) type. We have found that the transition 
is, in fact, of second-order (continuous). In the vicinity 
of the phase transition, the functional order parameter 
is given by the new scale-free solution that we described 
in Sec. IV Al This second-order transition is of a peculiar- 
nature, as it possesses critical exponents typically asso- 
ciated with first-order phase transitions. Consequently, 
numerical studies must be performed with care as finite- 
size effects can make this phase transition indistinguish- 
able from a first-order transition. 

We have found that in the vicinity of the transition 
the singular component of the free energy is F = 
where t — 7 — j c . The logarithmic correction is sufficient 
to make the first derivative dF/d7 continuous, therefore 
there is no associated "latent heat". Ordinarily, finite 
"latent heat" must be a continuous function of thermo- 
dynamic variables, which ensures that the phase transi- 
tion persists for finite T, at least up to the critical point. 
Conversely, when it is zero, it is plausible that the phase 
transition might disappear for arbitrarily small T as we 
claim. 

Throughout the paper we have attempted to keep the 
discussion as general as possible. All formulae derived 
using the replica method can be applied to a host of 
spin glass models defined on random hypergraphs. In 
the analysis presented in the paper it will usually involve 
the replacement of the cost function Ej(si, . . . , sk) by 
a suitable expression and performing disorder averages 
(. . .) j appropriately. Quantum analogues of dilute p-spin 
[13, if-NAE-SAT Bp, the Exact Cover [H, [H, Ha| , and 
the Vertex Cover [6(j] problems can be studied using this 
method. 

We have also devised a new method, of Quasi Monte 
Carlo variety, for the numerical determination of the 
functional order parameter. Since it significantly outper- 
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forms standard Monte Carlo, it can be used to improve 
the accuracy in the numerical studies of 1-step replica 
symmetry breaking, which so far required significant nu- 
merical effort f6ll | . 

For future work it is of interest to investigate the stabil- 
ity region of the replica symmetric solution on the plane 
(7, T) in a quantum regime corresponding to finite values 
of T. In the classical case, T = 0, the replica symmetric 
solution loses the stability at the point of the dynamic 
(replica-symmetry breaking) transition 7 = 7^. Beyond 
this point the energy landscape is characterized by a pro- 
liferation of an exponentially large (in N) number of deep 
local minima in the energy landscape, which traps clas- 
sical annealing algorithms. It is of interest to explore 
how this picture is modified for finite values of T. The 
structure of the free energy landscape will have implica- 
tions for the scaling of the minimum gap in the Quantum 
Adiabatic Algorithm. 

The effective classical Hamiltonian lj44j) may be used as 
a starting point for performing RSB analysis. Although 
it reflects static approximation, the disorder dependence 
is explicit. The replica-symmetric ansatz that we made 

I 



corresponds to the assumption the distributions of mag- 
netizations on different sites are not correlated. In the 
limit of small T relevant local minima correspond to in- 
teger values of magnetizations: m ~ ±1 and m ~ 0. 
Although local minima with intermediate values of m 
exist, our analysis indicates that corresponding free ener- 
gies correspond to excitations with energies much larger 
than the typical 0(T). In this limit, continuous magne- 
tizations may be replaced by discrete variables taking 3 
possible values. The third possibility (m = 0) makes the 
problem distinctly different from the classical if-SAT. 
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APPENDIX A: DETAILS OF REPLICA CALCULATIONS. 

1. The free energy functional. 

Substitution of the classical Hamiltonian of random if-SAT (fT4| into Eq. lfl7|) leads to the following expression for 
the disorder-averaged n-replica partition function: 

(Z n )= E e^M WW / J] cxp('- Cn ... 4A .E/ /3 dti?(< 1 W,...,^W;J n ... 4K ))\. (Al) 

[{«?(*)}] \ii<-<iK ^ a Jo 'I 

To perform the disorder average over c% x ,..i K , we expand the exponential and exploit the fact that averages (c? ■ ) = 
K\"f/N K ~ 1 are independent of p > 1. The resulting expression is 

(Z n ) = V e" JVWavg[{s?(t)}I+i: <"' c:[s?(t)] , (A2) 
[{»?(*)}] 

where we the term exp(— N7i avg [. . .]) results from disorder averaging. 

*WK(*)}] = U £ E^((Er dt ^(4(*)-^L(i))) ? ) • (A3) 

Like everywhere else, the disorder dependence of the cost function Ej({si}) is indicated with a subscript whenever 
J is a dummy variable and does not refer to a specific constraint (cf. Ji 1 ...i K in Eq. (|Aip ). In this particular case, 
an average over 2 K possible realizations of J is performed, which is indicated by the use of angle brackets (. . . )j. 
Similar to K-SAT, for many other combinatorial problems defined on random hypergraphs that are described with 
the present formalism, this notation is used to denote averaging over a "non-geometric" component of the disorder. 
For instance, for random if-XOR-SAT problem [52| the non-geometric component of disorder is a scalar J, taking 
values ±1 with probability of 1/2. However, for random Exact Cover [g] the disorder is purely geometric and the 
application of (. . . ) j produces no effect. 

As the next step we rewrite Eq. (| A3[) as a sum of terms with Ising-like interactions among spins (i.e. so that terms 
only involve products of spin variables). We observe that any function of K binary variables can be written as a 
polynomial in s\ . . . Sk ■ 

p ( a a \ - ptt«i}) «i a n K _ p(00...0) p(10...0) p(01...0) p(ll...l) 

{rai}e{0,l}* 
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In particular, for if -SAT, Eq. (|A4j) holds with Ey 1 '"™*' = 2 k-i J" 1 • ■ ■ ■ Using this transformation and inverting 
the orders of summation and integration we can rewrite Eq. (|A3|) as follows: 

P=l ai ...a p J {nrfJelO,!}?^ V=l 'J i!<— <i K 1=1 { r \n r i = l} 



- Havg =K\ 1 j2- z ^r- e / ^■■■dt p £ / n ^,..n^\ j_ E n n -w- 



(A5) 

From the form of Eq. I|A5|) , it is obvious that the partition function (|Al[) may be evaluated using mean field theory. 
The expression for Hn Vg [{s^ {t)}] can be written entirely in terms of Q ai ...a p (ti, . . . ,t p ) defined by Eq. (fTS)) . We may 
write H avg [{sf (t)}} = r ynl3V[{Qa 1 ...a p (ti, ■ ■ ■ where the newly introduced "quasipotential" is expressed in terms 
of spin correlation functions: 

1 00 1 f I p \ K 

V = ^E^i E E / d*i--d*p(n^i- -nrAJ)) HQ{a r \ nrt =iy({tr\nrt = l})- (A6) 

P P=1 P ' {n rl },...,{n rK }e{0,l}Pa 1 ...a p J V=l ' J 1=1 

The only other term that contributes to the replica free energy, is of an entropic nature, as it appears from 
the summation over all spin configurations subject to the constraints lfl8|) on correlation functions. This term 
will be referred to as "quasientropy" and denoted S. In order to compute the quasientropy, we introduce auxil- 
iary variables A ai ... (ti, . . . , t p ). Each constraint of the form Q = -^J2i x i ls enforced by inserting the integral 

/+^ d\e~ NXQ fli eAxi - We ma y write 

S = iV^3 ln / VXe ~ N12pl2ai ■" fitl '" div Aai - ■»^ tl -- t ")«°i- ■^ {tu -' tp) Z?l{\ ai ... ap (t 1 , . . .,t v )}}, (A7) 

where T>X denotes a functional integration over the set of X ai ...a p (ti, ■ ■ ■ , tp), and Z\ is the "single-site" partition 
function 



2 l[{A}] = ^2 eX P(E E / dt l--- dt P X ai...a p (t 1 ,...,tp)s ai (ti)...S ap (t p )+^2lC[s a {t)} 
\{sa(t)}] \=la 1 ...a p - ) a 



(A8) 



The partition function may be written in the form of the functional integral over both Q ai ...a p (ti, . . . ,t p ) and 
X ai ...a p (ti, ■ ■ ■ ,t p ) as seen in Eq. (fl9|) . This functional integral is in the limit N — > oo dominated by its saddle- 
point value corresponding to the extremum of the free energy functional 

F[{Qh {A}] = lV[{Q}} - S[{Q}, {A}]. (A9) 

Normalization constants have been chosen so that the real free energy F = (—l/Nnf3)ln(Z n ) corresponds to the 
extremal value of IIA9I). 



2. Static Approximation 

The saddle-point solution may be obtained by varying the free energy functional •? r [{Q},{A}] with respect to 
its time dependent arguments; and self-consistency equations for saddle-point values of Q ai ... ap (ti, . . . ,t p ) and 
A ai ...a p (ii) • • • ,tp) may be written down. Unfortunately, the resulting infinite number of equations cannot be re- 
duced to the closed-form expression for the free energy. To have a closed-form solution is important since an analytic 
continuation to n — > must be performed eventually. 

We resort to an approximation of variational type, which restricts the functional integration in l|19p to only time- 
independent Q ai ...a p and A ai ... Qp . In other words, the functional integral is replaced by a regular, albeit infinite- 
dimensional integral. Under this static ansatz, the saddle-point condition for the integral over X ai ...a p becomes 

where the static single-site partition function is rewritten in the following form: 

Zl ({A 0l ... 0p }) = eX p(E E Aa i- ■ a P / dtl " ' dt P S «l(*0 -S a P (*p) +^2 IC i S a(t)} 

K>a(t)}] \=l ai ...a p J a 



(All) 
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Stationarity condition (|A10|1 shows that the saddle-point values of the static order parameter Q ai ...a p are the time- 
averaged correlation functions 

Q ai ...a P = l^r [ dt 1 ---dt p s ai (t 1 )---s ap {t p )) . (A12) 



Here (■■■)z 1 represents an average taken with respect to the Gibbs distribution corresponding to the single-site 
partition function l|All|) : 

p[{s a (t)}] = i-e^-i"^ 1 •■-p/ d *i- d Ap'.i(*0"-.,(*p) i (A13) 



Once time-dependence is ignored, all expressions may be rewritten in terms of the mean magnetizations m a — 
4 Jq dt s a (t) . In particular, we introduce the marginal p({m a }) of the Gibbs distribution 1|A13[) 

P(K})= E P[{s a (t)}}l[5(m a -U P dts a (t)). (A14) 

[{-.(*)}] - V ' Jo J 

Note that for finite L the delta-function is the Kronecker's delta; the normalization has been chosen so that 
/ {dm a } p({m a }) = 1. We replace sum over paths [s(t)] by integrals over magnetizations. Eq. (|All[) may be rewritten 
as follows: 

Z x = / {dm.Je^/^.-p^-^" 1 ""-''^" 1 ''"-', (A15) 
where the factor e~P u o( m ) comes from the summation over paths of average magnetization m 



e -/3u (r. 



= Ve^Wm-i/" dt Sa (t)). (A16) 
[.(*)] V PJo 1 

For a finite value of L it may be written as a sum over kinks of s(t). In the limit L — > oo the series may be expressed 
in the analytical form given by Eq. I|30p . 

We substitute Q ai ...a p = J {dm a } p({m a })m ai . . .m ap into Eq. I|A6|) and gather all terms in an infinite sum to form 
an exponential. The resulting expression for the quasipotential reads 

Vb(W)] = ^J {dmW}p({™«}) • • -p({mW}) (l - e-"£. , (A17) 

where Ej(m\, . . . , rriK ) = X){ ni }e{o i} K Ej^rn^ 1 • • • m 7 ^ is the generalization of Ej{s\,...,sk) to continuous 
magnetizations. 

The expression for the quasientropy, expressed in terms of p({m a }), is 

S[p({m a })} = -J {dm a } p({m a })(\np({m a }) + pJ2 u o(m a )) ■ (A18) 



3. Replica-symmetric ansatz and the functional order parameter 

Now that the free energy functional corresponding to the static approximation has been written down, we must 
perform analytical continuation in n and take the limit n — > 0. To accomplish that we employ the method developed 
by R. Monasson for systems with binary variables [20l . l2ll . l63l | and extend it to continuous magnetizations. 

Replica symmetric ansatz assumes that the saddle point of the free energy functional, the distribution p({m a }), is 
a symmetric function of its arguments. We make the following ansatz for p({m a }): 

r e -Ph(m a ) 

p({m a }) = [dh(m)] P[h(m)] ]J fdme _ ph[m) . (A19) 

Here P[h(m)] is the probability distribution of functions h(m) — the "effective fields" that define Boltzmann dis- 
tribution of magnetizations. It is evident that p({m a }) is symmetric and normalized to unity, provided that 
f[dh(m)]P[h(m)] = l, 
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Substituting (| A19|) into (|A17|) we obtain for the quasipotential 

V= -J [{d^(m)}f =1 ] P[h 1 (m)]---P[h K (m)}H {dmW}f = ° " ' (l - C -^A^\-^) 

(A20) 

Product over n replicas appearing inside the functional integral may be replaced by its argument raised to n-th power. 
Upon taking the limit n^Owe use x n w 1 + n In x to obtain the expression given by Eq. f25|) . 
For evaluating the quasientropy l|A18(l we use the following trick described in Ref. [63 |: 

d f 

{dm a }p({m a }) lnp({m Q }) = lim — / {dm a } (p({m a })) r , (A21) 



♦ l dr 

where r-th power of p({m a }) for a positive integer r can be written in the form 

f ^ p -/3(fti(m )+-+/i r .(m a )) 

(P({m a })) r = J [d hl (m) ■ ■ ■ 6hr(m)]PMm)] ■ • • P[h r (m)] ]J ^ AmG _ [jhl[m)) ; ; ; ^ dwe _^ M) • (A22) 

In the limit n — » products over n replicas may be replaced by the logarithm so that the r.h.s of Eq. (| A22|) splits 
into several additive terms. It may be further further simplified to 

{dm a } (p({m a })) r KnJ [dh(m)] (F H [/i(m)] - P[h(m)fj In / dme- 0h ^ m) . (A23) 

We write P^* T ^[h(m)] to denote r-th functional convolution of P[h(m)] with itself: 

P l * r] [h(m)} = ( P*P*...* P) [h{m% (A24) 

r times 

where the convolution of two functionals A[h(m) and B[h(m)) is defined as follows: 

(A * B)[h(m)\ = J [du(m)} A[u(m)]B[h(m) - u(m)]. (A25) 
The convolution P^* r \m) can be written as an inverse Fourier transform of (P[w(m)]) , where 

P[u(m)] = [ [dh(m)] cT^ dm«(m}h(m) p[h( m )} . (A26) 



Once Eq. (|A23[1 is expressed in termos of (-P[w(m)]) r , it is straightforward to perform an analytical continuation in r. 
We now turn to the remaining term in (| A18[) : 

{dm Q }p({m a })]>> (m a ) =n J [dh(m)] P[h(m)] J ^£^ggff 

= „ / [dMm)] ^^^;!^! 7 ^ / [<M»0] e^-^PKm)]. (A27) 

We substitute 

„— /3/i(m) i x /• 

rj , ^7 ,, = - ~ 7^ In / dm' e~^ m ) ( A28) 

/ dm' e-Wro ') /3 <5/i(m) 7 v y 

into (|A27P and integrate by parts. 

Putting together (|A23P and (|A27j) . expressed entirely in terms of P[ui(m)], we obtain the expression lf28|) for the 
quasientropy. 
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4. Classical limit (T = 0) and longitudinal stability 

In the limit r = the expression lj30|) reduces to 

e-^™) =5(m-l)+S{m + l). (A29) 

Since all integrations over magnetizations necessarily involve a factor of e~P u °( m \ they can be replaced by sums over 
to = ±1. Specifically, we will write 

e -0h(m) _ e -Ph+is(m - 1) + e-' 3 ' l - 1 ( 5(TO + 1). (A30) 

The free energy functional JF[P[/i(to)]] can be expressed in terms of a joint probability distribution of h + i and h-\. 
It is convenient to make a change of variables h±i — g^fh. Expressed in terms of the probability distribution P2(h, g), 
the free energy functional has the form 

K 



T=l[{dh e dg e } f[P2(h e ,ge) x In + e -^'+ h ^ - ln^ . 



l=\ v i {s e } 

- y dhdg In (e-^ + e'^) J ^ e^ h+i ^P 2 (u, v) (l - In i/)) , (A31) 
where P2(uj,v) is the Fourier transform of P2(h,g). After some simplifications it becomes 

JF=7 / {d/i,} JJ P (^) x fe ln ( 2cosh ^) -biJ2 e ' l3Bj{Sl '"'' SK)+PJ:eheSe ') 
e=i ^ e {s e } ' 

-ljdh\h\J^ e^ h P(u) (l - lnP(u;)) + AF[R(v)], (A32) 

where P(/i) = / dg P-2{h, g), the corresponding Fourier transform is P{oj), and the last term is a functional of R{v) = 
ft (0,1/): 

Ajf [£(„)] = ljdggj^- e^R(v) (l - In %)) . (A33) 

This contribution is zero, provided that the inverse Fourier transform of lnR(v) is non-negative everywhere except 
the origin. This will happen for the saddle-point of T\ incidentally, it ensures that R{g) = fdhP2(h,g) 0, as should 
be expected for a probability distribution. Without the last term, Eq. (|A32|1 coincides with the free energy functional 
obtained in Ref. 12111 - 



We can express the condition for the longitudinal stability (that is, the stability within the replica-symmetric sector) 
of P(h) that makes l|A32|) stationary as the requirement that eigenvalues of the Hessian matrix d 2 T '/ ' dQ ai ... ap dQb 1 ...b q 
associated with replica-symmetric eigenvectors be positive. In the classical limit (T = 0) the number of independent 
order parameters {Q{k r }} is greatly reduced because Q{k r } = Qp for p = X] r &2r+i- As a consequence, the free energy 
may be expressed entirely in terms of the parameters OfciOO..., which means that we may assume that a\,...,a p are 
all different (as well as b\, . . . , b q ). 

The condition that the Hessian is semi-positive definite can be written as a condition that for an arbitrary {u p } p *L 1 

Sp.q U P U q Sa 1 ...q p J2 bl ...b q g W/dQg, . . .„„ 8Q bl . . . bq 

V u 2 T 1 ^ 6 > 

Note that retaining the denominator in l|A34[) is mandatory. This is peculiar to the replica theory, as in the limit 
n — > the denominator may become negative. We can choose {u p } in the following general form: 

dh mt) u{hy (A35) 
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Substituting this expression along with the degeneracy factor ^2 ai a 1 = (™) into (|A34p . we obtain 

/ ^^WSWs? C)»? > »■ < A36 > 

The denominator includes only terms with even p, provided that u(—h) = u(h), and, conversely, only terms with odd 
p when u(—h) = —u{h). In the limit n — > the binomial (™) « (— i.e. the denominator is negative for even 
u(h) and positive for odd u(h). 

Therefore, the stability condition is that the free energy functional T[P{h)] is a local maximum with respect to 
even perturbations and a local minimum with respect to odd perturbations. In the cases where the symmetry of the 
problem at hand permits one to restrict the space of P(h) to even functions, as if -SAT does, one can safely maximize 
the free energy within this restricted subspace. 

In a quantum case, the stability condition may not be expressed in such a simple form. Fortunately, for if-SAT in 
a quantum limit, the stationarity condition determines saddle point uniquely, obviating the need for analysis of the 
longitudinal stability. 



APPENDIX B: CAVITY METHOD AND THE BETHE-PEIERLS APPROXIMATION 

The drawback of the method of replicas is the lack of intuitiveness. However, all results obtained using replica 
theory can alternatively be derived using the so-called cavity method that does not rely on analytic continuation in n. 
In fact, many physical properties are more easily derived using the cavity method. It was recentlyused by M. Mezard 
and R. Zecchina to obtain the 1-step RSB solution to classical zero-temperature random if -SAT I22j.[2al . An excellent 
description of the method is given in Ref. [24]. In contrast to the replica method, cavity equations [67] are written for 
a specific realization of disorder; the disorder averaging is performed as a last step. At the level of replica symmetry, 
cavity equations are identical to belief propagation (BP) equations as was demonstrated in [23 |. 

We generalize BP equations to include quantum degrees of freedom. Although we are limited to static approx- 
imation, this excercise allows us to establish the physical meaning of the former. By examining the form of the 
self-consistency equations f36|) for the order parameter, we conjectured a particular form of the effective Hamiltonian 
given by Eq. (|31a[l . The static ansatz is prominent in that the effect of individual spins in the effective Hamito- 
nian is only via their magnetizations, viz. imaginary-time averages of the spin. We corroborate the conjecture by 
demonstrating that we rediscover the results obtained in Sec. Hill using replica-symmetric static ansatz. 

We introduce the Gibbs distribution p(mi, . . . , mjv) associated with partition function I131I) 



p({mi}) = z(lJ}) e ^ Cff({ms};{ ' 7}) - (B1) 

The probability distribution of magnetization of a particular spin Pi(m) may be written as a marginal of p({nii}), 
where all magnetizations other than rrii have been integrated over: 



Pi(m) = j {&m v } v ^ip({m v \mi = m}). (B2) 

We write p({m,i> jm, = m}) to indicate that the z-th argument of p(m\, . . . , mjsr) should be set to m. 

We work with the ensemble of random hypergraphs, where each hyperedge may appear with probability Klj/N K . 
Randomly chosen vertex i will have on average if 7 hyperedges incident to it. The degree d{ of the vertex is Poisson- 
distributed with mean if 7 (see Eq. ©). We temporarily relabel vertices as follows (see Fig. [22]) : The randomly 
chosen vertex shall be referred to as the cavity vertex and shall be labeled by number 0. Its immediate neighbors shall 
be labeled by two integers: k = 1, 2, . . . , d = dj, to refer to one of the incident hyperedge, and I = 2, . . . , if , to label 
other vertices connected to the cavity vertex via this hyperedge. The particular order of variables for each hyperedge 
or the ordering of hyperedges is unimportant due to symmetry; however, we must remember to reorder components of 
Ji 1 ...i K appropriately. We shall write Jk to denote disorder variables associated with fc-th (k — l,...,d) hyperedge; 
the reordering of components of Jk is such that the combination of variables with sq = jj. 1 ^ and Ski = jjf for 
I = 2, . . . , if , is penalized. 

The following is the description of the Onsager's cavity method. We use to denote the partition function of a 
system with N — 1 spins obtained by removal of a cavity spin. The partition function Z of the original system can be 
related to Z^ via 

Z = Z ( - a) [ dm{dm ke }e~P uo{m) -P^=^ {m < mk2 --- mkK - Jk) . (B3) 
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FIG. 22: Demonstration of the cavity method. The white circle is a cavity spin, and the white triangles are the hyperedges 
incident to it. The black circles are the neighbors of the cavity spin. p^}{m) denotes the probability distribution of magne- 
tizations of neighbors realized if the cavity vertex and hyperedges incident to it are removed. The probability distribution of 
magnetizations of the cavity vertex po(m) is completely determined by {p^e ( m )} as we indicate with gray arrows. 



Since the cavity vertex has been chosen at random, the structure remaining after the removal of the vertex and its d 
incident hyperedges is a random graph of the same average connectivity AT 7 with N—l vertices. Moreover, locations 
of the vertices labeled by k and £ are random and uncorrelated in the absence of a cavity spin. It is known from 
the theory of random graphs that randomly chosen vertices are either unconnected to each other or separated by at 
least O (log AT) edges. As a result, distributions Py>{ m ) realized in absence of cavity spin are uncorrelated for different 
k,£, unless a long-range order is present in the system. The absence of this long-range order is equivalent to the 
assumption of replica symmetry. 

Using (|B3|) we can relate the distribution po(m) of magnetization of the cavity spin to a set of pj^(m): 

po( TO ) = const x I {dm M }e- ,J "» (m) ^ E '=^ (ra '"" ! -- ra » c;J ' ) Y[Pkt( m kt) 

k,t 

d . K 

= const xe^ Uo(m) JJ / {dm*} e -f>E(™,m a ,..., m K -,J k ) JJpg 5 ("**)• (B4) 

k=l 1=2 

We perform the change of variables hi(m) = const — A lnj»j(m) to rewrite (|B4[) in the form 

d 

ho(m) — up(m) + ttfc(m) + const, (B5a) 

fe=i 

Ufc(m) = u Jk (m; [h^(m), . . . , h^K(m)]), (B5b) 

where the function uj(m; [{hi(m)}]) has been defined in l[34|) . 

To write a closed system of equations we employ the fact that the probability distribution of a disorder in an A^-spin 
system factorizes. The number d of hyperdedges incident to cavity vertex and vectors {Jk}k=i are n °t correlated 
with the disorder distribution of the remaining (N — l)-spin system. To each disorder realization there corresponds 
an effective field ho(m) associated with the cavity vertex. Accordingly, we may speak of the probability distribution 
of effective fields P^ N ^[h(m)]. Absent long range order, effective fields h^(m) associated with neighbor vertices are 
independent and drawn from the distribution P^ N ^) [h(m)] for a system with A^ — 1 spins. With the aid of auxiliary 
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distribution Q N - 1 [u(m)}, we can relate P (JV) [/i(ra)] to P^ N -^[h( 



m 



K 

QW- 1 \u(m))= {dh e (m)} Y[P<- N -^[h(m)} x (5[u(m) - Uj (m; [h 2 (m), . . . , h K {m)])\) j, (B6a) 

J 1=2 

oo „ d i- d 

pW[h(m)}=J2fd(K7) {dh k (m)} [] P^[h k (m)} x 6 h(m) - u (m) - 5> fc (m) 

<i=0 fc=l L fc=l 



(B6b) 

In the asymptotic limit TV — > oo we may replace p W ffe (m)] and P^ -1 ) [/i(m)] by the same limiting distribution 
P[/j(m)] recovering self-consistency equations (|36a,|) . (|36bl|) . 

The free energy for the concrete realization of disorder can be obtained using Bethe-Peierls approximation [62]. 
This approximation is exact below the percolation threshold and should be asymptotically correct everywhere in the 
replica-symmetric phase. Bethe free energy takes the form of a functional that depends on the probability distribution 
{Pi(m)}fL 1 as well as the joint probability distributions Pi^...i K ( m i, ■ • ■ i m K) associated with hyperedges. 

For the remainder of this section we shall write Yl(i 1 ...i K ) ■ • ■ to indicate sum over all hyperedges, alternatively 
written as J2i 1 < - <i K c ii-..i K x • • • . For example, the degree of vertex i can be written as d± = £V iK <. &it,i- The 
Bethe free energy functional associated with static Hamiltonian (|31a|) can be written as follows: 

•pBcthc[{Pi("l)}, {Pu...i K ( m l) ■ ■ ■ i m K)}] = 

J {dm^}jjj 1 ..., y {{me})(^lnpi 1 ...j K ({m£})+E({me}]Ji 1 ...i K j)-y^ y (d i -l) J dmpi(m)(^kipi(m)-uo(mfj. 

(i 1 ...i K ) i 

(B7) 

This free energy is to be minimized, subject to the constraint 

Pi e (m) = J {&m£i}e^tPi 1 ...i K {{rrii'\rni = m}), (B8) 

as well as normalization conditions J '{dmj pi l ...i K { m ii ■ ■ ■ i m K) — J dmpi(m) — 1. Introducing K X M Lagrange 
multipliers f associated with constraint (|B8[) and extremizing the Bethe free energy, we obtain the following 
equations to be solved self-consistently: 

K 

Hi(m) = u (m)+ ^ ^ dii,i u i7...iK ( m ) + const, (B9b) 
(ii. ..»«•) «=i 

^L. (m) = H u (m) - uV Ak (m). (B9c) 

The right-hand-side (r.h.s.) of Eq. (|B9b|) is the sum of uo(m) and contribution from all di hyperedges incident to 
vertex i; r.h.s. of Eq. (|B9c|i is equivalent to the sum of Uo(m) and di t — 1 hyperedges incident to vertex it with the 
exception of the hyperedge (ix, . . . ,ik) (see Fig.[23]for the illustration). 

These self-consistency equations obtained from Bethe-Peierls free energy functional are known as Belief Propagation 
(BP) equations in statistics [62], or as Thouless- Anderson-Palmer (TAP) equations in spin glass theory jjj]. The BP 

algorithm updates local "beliefs" {uf^ A (m)}, {hf^ A } until convergence. These beliefs describe distributions Pi(m), 
Pii...i K (TOi, . . .,m K )\ 

Pl (m) = const xe-^-f" 1 ', (BlOa) 
Pil ,„ iK (mi, ...,m K ) = const xe -^>.-,>»ir!J (1 ..., lf )-^L<... jjf W i (BlOb) 



In the regime of the stability of the replica symmmetry, solutions to BP equations (|F39[) should be unique with high 
probability. 

The disorder average of the Bethe free energy can be derived in terms of P[h(m)] and Q[u(m)] — the histograms 

°f {^ f.iy } and {uf} That P[h(m)] and Q[u(m)] must satisfy self-consistency equations ([361 follows directly 

from i|B9all and (|B9c|l . The histogram of Hi(m) is also P[h(m)]. 
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FIG. 23: Illustration of the BP equations for (K = 3). The white triangle is a hyperedge ■ With this hyperedge we 

associate 3 functions (beliefs) u ( i 1 i \ 2i;j (m), i4fi 2 i 3 ( m )j and Wi^ 2 j 3 (wi). Beliefs {^i*i 2 , 3 } correspond to sum over all hyperedges 
incident to it. other than {iiiiiz). Grey arrows indicate that i4*]^(m) is determined by /4ii 2 i 3 ( m ) an< ^ '^i^is^' 71 )' Similarly 



M^ 2i3 (m) is expressed in terms of h\^ i2i3 (m) and ^^213 ( m ) : an d so on 



(i) 



,(3) 



We shall require another distribution: the histogram of {Hi(m)} but with each vertex weighted in a proportion to 
its degree dj. With Pd[H(m)] denoting histograms of Hi(m) for vertices of a degree d, the desired distribution is 



P'[H{m)] = ^Y, d fd(Ki)P d [H{m)\. 



Due to the fact that the distribution of degrees is Poisson, Eq. (|B11[) may be rewritten as follows: 

P'[H(m)] = J [du{m)]Q[u(m)]P[H(m) -u(m)]. 
Substituting IjBlOaj) and (|B10b|l into (|B 7|) and performing the disorder average, we obtain 



(Ml) 



(B12) 



(F 



Bcthc/ 



J[{dhi(m)}} f[PMm)} x (in J {dm e }< 



-0E({m e };J)-/3 T,f=i M«w) 



+ ^-J [dh(m)du(m)] P[h(m)]Q[u(m)] In J dme~' 3h{m) -' 3u{m) -^J [dh(m)] P[h[m)] In J dme^W. 

(B13) 

This coincides with the expression obtained using the replica method [see Eqs. lf25|) and 1(38])]. 



APPENDIX C: QUASI MONTE CARLO IMPLEMENTATION 

The Quasi Monte Carlo (QMC) method of evaluating integrals replaces the random sequences of standard Monte 
Carlo (MC) algorithms with deterministic minimum discrepancy sequences (53|. For example, a two-dimensional 
integral of a function f(x, y) on [0; l] 2 is approximated by 



N-l 



[Oil]' 



dxdyfix^n^-jY,* 



(CI) 



k=l 



where {ik, jk}k=i i s the minimum discrepancy sequence (since it is a finite, we will call it a minimum discrepancy 
set). We use Sobol sequences and choose N to be an integer power of 2 for best results. Error estimate for 
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the two-dimensional integral (|C 1[) is 0(\ogN/N) for continuous functions f(x,y) and 0(1/7V 2 / 3 ) for discontinuous 
functions f{x,y). This compares favorably to the expected error of 0(1/N 1 / 2 ) in standard Monte Carlo. 

We adapt QMC to integrals involving univariate probability distributions. The probability distributions will be 
represented internally as a finite-size sample. In contrast to the standard MC, we will ensure that these samples are 
as uniform as possible. With each univariate distribution p(x) we associate a function X(p) defined on the interval 
[0; 1] and satisfying the condition 



X \J dx'p(x')j = x. (C2) 

Internally, it will be represented by a set of {^A;}^^ 1 , where Xk — X(k/N). 
For an arbitrary function X(x), its expectation value may be approximated by 

dxp(x)\(x)= 1 ^- Y J2 X ( X ^ ( C3 ) 

The flow of the computation shall consist of a sequence of transformations of probability distributions. The elementary 
operation is finding the distribution of a variable z = f(x,y), given distributions of variables x and y 

p(z) = dxdyp(x)p{y)S(z - f{x,y)). (C4) 



That is, we need to find a uniform sample {Zk} of a distribution p(z) from uniform samples {Xk} and {Yfc} of 
distributions p(x) and p(y). What the appropriate sample should be, can be assessed indirectly by considering the 
expectation value of an arbitrary function X{z): 

N-l . 

— £ \{Z k ) = / dzp(z)X(z) 

i i « 



N 

fe=l 



'[Oil]: 

Estimating the integral over [0; l] 2 in (jC5j) and using 1|C1|) . we may write 



dxdy p(x)p(y)X(f(x,y)) 

d Pl dp 2 X(f(X( Pl ),Y(p 2 ))). (C5) 



W^J E = tv^I E A (/(^ fc , Y jk )), (C6) 

k k 

where {ik, jk}^=i ls the Sobol set. The choice of the sample {Zk} satisfying <|C6|) for any X(z) is unique. Algorith- 
mically, it is computed as follows: 

1. For k = 1, . . . ,N — 1 evaluate Zk = f(Xi k ,Yj k ), where {ik,jk}k=i i s the Sobol set. 

2. Sort the resulting vector {Zk}k=\ m increasing order. 

The last step is to ensure that Zk < Zk+i, which is required by definition (|C2|) . 



1. Application to T = classical 3-SAT 

Let us briefly describe how this idea can be applied to solving self-consistency equations. For K — 3, Eq. (|65a,|) 
already has the form of l|C4p . 



Q{u) = J dh 2 dh 3 P(h 2 )P(h 3 )5(u + J 1 mm(l 7 {J 2 h 2 ) + ,(J 3 h 3 ) + )), (C7) 
enabling one to compute {uk} from {hk}. Eq. <|65b|) can be rewritten as follows: 

oo 

P(h)=J2fd(^)Pd(h), (C8a) 



d=0 



p d (h) - f {duk} n Q(u k ) xs(h~Y, u k)- ( csb ) 

J fe=l ^ fe ' 
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A set of {Pd(h)} is computed using the following recursive definition having the desired form of l|C4j) : 

P d (h)= { Ah'AuP d - X {h')Q{u)5{h-h' -u), (C9) 



together with the condition Po(h) = S(h). The distribution Po(h) is represented by a vector of N— 1 zeros. Computing 
a sample of P(h) from a set of samples of Pd(h) via (|C8aj) means that we have to select N — 1 values from the larger 
set of (dmax + 1) x (iV — 1) values that represent distributions Po(h) through -Pd max (^)- 

The elegant way to accomplish it is the following. We formally introduce the function H(t,p) defined on [0;1] 2 . 
For any fixed value of t, viewed as a function of one argument p, H(t,p) represents the distribution Pd(t)(h): 



h(J dh'P m (h')\ =p, (CIO) 



and d(t) is a step-wise function of t such that 

d(t) d(t) + l 

^/ fe (3 7 )^t< J2 mw- ( cn ) 

The expectation value of an arbitrary function \(h) can be written as 

[ dhP(h)X(h) = [ dtd P X(H(t,p)). (C12) 

J J[0;1] 2 

Applying IjCip to the integral, we construct the sample for P{h) from {h^* k ^}, where {h^} represents the sample 
for P d (h). 

Memory requirements for each iteration step can be kept at O(N), whereas the time complexity is 0(N log 2 N), 
which is the product of d max = O (log AT) and the 0(N log N) complexity of sorting. 

The procedure described above is trivially extended to K ^ 4 and to finite temperatures T > 0. The extension 
to finite temperatures merely changes the form of Uj(h 2 , ■ ■ ■ , h K ). Distribution Q(u) may be computed using either 
a single (K — l)-dimensional integral, or as a sequence of K — 2 two-dimensional integrals. The latter approach is 
possible because for any K ^ 4, the function uj Q12, • ■ ■ , Hk) may be written in terms of compositions of functions of 
two variables. 

Instead of iterating self-consistency equations for a fixed value of 7, we achieve accelerated convergence by specifying 
the desired width A of the distribution P(h) and adjusting the value of 7 at each iteration step to satisfy this constraint. 
As a result, we observe exponentially fast convergence and avoid the effects of the critical slowing down in the vicinity 
of the phase transition. 



2. Application to quantum K-SAT 



The case of quantum K-SAT in the limit T -C 1 is slightly more involved. The order parameter is the joint 
probability distribution (j.p.d.) P(h,h). For Quantum Model B, it is possible to parameterize this j.p.d. by a 
univariate distribution P+(h) and apply the method described previously. For Quantum Model AB, however, no such 
parametrization is possible. 

The workaround is to work with univariate distributions Q(u) and R(t]) exclusively. It is straightforward to compute 
Q(u) from R{q) by evaluating the {K — l)-dimensional integral (|6Tj) that can be further reduced to the sequence of 
two-dimensional integrals for K ^ 4, The non-trivial part is the computation R{rj) from Q(u), which we describe 
below. 

It is easy to see that Q{u) is symmetric [i.e. Q(u) = Q{— u)] and that |u| < 1. Let £i/ 2 denote the probability that 
M>l/2: 

£1/2=/ duQ{u). (C13) 

J\u\>l/2 

We introduce the reduced distributions Q<(u) and Q>(u) defined on intervals [0;l/2] and [— 1; —1/2] respectively: 

<9<W = r4— Q(u)6{u)6(\-u), (C14a) 
1-51/2 V2 ' 

Q > (u) = -^Q{u)e{-u)e(-\-u), (C14b) 
£1/2 v 2 / 
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I 1 for x > 

where 6{x) is the Heaviside function 6(x) — i '. The factors 2/(1 — £1/2), 2/£ 1 / 2 ensure that Q < (u), 

10 for x < 0. 

Q> (u) are normalized to unity. 
We also define 

W) = f>(^(l-£i/ 2 )) f{du k } YlQAu^xsfh-J^uX (C15) 

d=0 ^ k=l ^ k ' 

as well as the sequence {P/-(ft)}: the sequence of successive convolutions of Po(h) with Q > (u). It is computed via the 
recurrence relation 

P k (h)=JduQ > (u)P k _ 1 (h-u). (C16) 

The distribution Pk(h) gives the contribution from vertices that have arbitrary number of incident hyperedges with 
b6 [0; 1/2] and precisely k hyperedges with u € [—1; —1/2]. In view of a relation (|55[) between u and u, to each h in 
the distribution Pk(h) there corresponds h — k + h. 

Including contributions from mirror image regions u S [— 1/2; 0] and u G [1/2; 1], the distribution P(h,h) may be 
written in terms of {-P/-(/i)}: 

p(h,h)= fk+^Yti/*)/*- (^^1/2) / d^_A + (ft+)A_(M^-'»+ + M^-*M--*--'»+-M. 

(C17) 

where /fc(a) denotes the Poisson distribution with mean a as usual. It follows that the distribution R{rj) given by 
Eq. f97j) may be written in the following general form: 

R(rj) = J df + dt-d/i + dL P^it/i+lPitLjpi-Ji (i) - A( i+ ), k (L)Ci+, L)) , (C18) 

where fk+,k- C&+> ^-) = ? 7(^+ — /i-, T — fc + — fc_ — h + — hJ) and we have defined a step-wise function chosen to 
satisfy 

fc(t) „ fe(t)+i 

E^(^V»)<*< E *(^v0- (C19) 

r=0 r=0 

Once i?(ry) has been expressed as a 4-dimensional integral, values {^fe}^ 1 may be sampled using Sobol sets. Memory 
and time requirements of this procedure remain 0(N) and Q(N log 2 N), respectively. 
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